program mkdecofig; {Cornell July 99...}
uses common,d04_header,ppot,psdict,inp,outp,tools,decor;
Const
   tbtyps  = 20; {number of bond types}
   abond  = 20; {atom-bond types}
Type
   v4	= array[1..4] of myreal;
Var
   ea				: array[1..110] of myreal;
   nea				: array[1..110] of integer;
   i,j,k,l,o,oind,puckmode,ecol	: integer;
   txt,txp,txyz,txfig,txinp	: text;
   over,max_size,min_size,lwidth,bview,
   x1,y1,x2,y2,cparxyz,
   c1,c2,c3,c4,rx,ry,rx1,ry1,rx2,ry2,
   cc,det,epsw,cmp,zlim,zpuck:myreal;
   u,v,vv,w,r0,r1,r2,psw:v2;
   show_z_alab,show_atom_elab,show_z_alab_rel,show_uc,show_page,
   shift_view_t,shift_view_m,
   show_z,show_bonds,show_atoms,show_orb_labels,pboxstroke,
   show_atom_labels,show_atom_node_labels,show_occupancy,show_tbf,
   show_nu,show_objects,inv_int,show_node_num,show_node_mon,show_mon,
   show_node_flag,
   show_node_mon_mark,show_perp,show_tiling,show_tiling_perp,eofile,
   show_tbonds,show_atom_bonds,autoadv,autodouble,nobound,newpage,pse,
   title,asprb:boolean;
   addobj,show_xyz:integer;
   nadv,npsdefs,goback:integer;
   msg,ticoname,fxyz,tname,psout,caption : mystring;

   pswin,pswin_u,psclip	: v4;
   rec,zrec0,rec0,rece,olx:v2;
   nug : array[1..tbtyps] of i2;
   rho,atlw,atlwx:array[0..z_ele] of myreal;
   otlw:array[1..t_orb] of myreal;
   ots:array[1..t_orb] of boolean;
   ogrs:array[1..t_orb] of mystring;
   agrs,agrsx,atom_plot_str:array[0..z_ele] of mystring;
   axxx:mystring;
   ats,atsx:array[0..z_ele] of boolean;
   tbondthick,athick:myreal;
   dbg:integer;
   fsize,fsize0,zmin,zmax,z1,z2,rho_min,rho_max,cx,cy,aspr:myreal;
   dzmx,dzmi,demax,demin:myreal;
   ozmx,ozmi:integer;
   e1,e2,j1,j2:integer;
   u3,w3:v3;
   bw:array[-tfla..tfla,-tfla..tfla] of myreal;

   ch:char;
   tinp,loops,ii,ll,llp,lll,nfig,advm:integer;
   beta,ene:array[1..tloops] of myreal;
   sinfo:array[1..tloops] of mystring;
   tb:array[-tfla..tfla,-tfla..tfla] of boolean;
   obj_b:array[0..2,1..t_obj] of boolean;
   obj_lw:array[1..t_obj] of myreal;
   obj_fils:array[1..t_obj] of mystring;
   obj_start:array[1..t_obj] of integer;
   obj_fr:array[1..t_obj] of integer;
   x2y:boolean;
   rax:array[1..nsites] of v3;
   sene,dene:array[1..nsites] of myreal;
   sgs: array[1..nsites] of integer;
   ratp:array[1..nsites] of integer;
   nxyz,nats:integer;
   agrsb:array[1..abond] of mystring;
   nab:integer;
   ble:array[1..abond] of v2;
   chid:array[1..abond,0..z_ele,0..z_ele] of boolean;
   oid:array[1..abond,0..t_orb,0..t_orb] of boolean;
   nrows,ncols,irow,icol,nrc,show0,show1:integer;
   tbfr:v2;
   advance:array[1..tloops] of boolean;
   first:boolean;
   tnodes1,tnodes2,atoms1,atoms2,dbl,pages,tpages,m,n:integer;
   bxyz:array[1..nsites] of integer;
   zzr:array[1..nsites] of myreal;
   nnup:integer;
   nulist:array[1..5] of integer;
   n_ocr:array[1..tdim,1..at_of_node] of myreal;
   show_at_oc:array[0..z_ele] of boolean;
   oldraw:array[1..tdim,1..t_orb] of boolean;
   oldrawn:array[1..t_orb] of boolean;
   show_arrow:boolean;
   n_dz	: array[1..tdim,1..at_of_node] of myreal; {node.atom.z}
   n_dxy: array[1..tdim,1..at_of_node] of v2; {node.atom.z}
   pp : array[1..nsites] of i2; {for draw_xyz}
   dsca:myreal;
   xx:real;

procedure stacking_faults;
var
   j,l,o,k,io : integer;
   nsf	   : array[-tfla..tfla] of myreal;
begin
   if(multi>1)then io:=1 else io:=0;
   for o:=-tfla to tfla do nsf[o]:=0;
   for l:=1 to tnodes do begin
      for o:=-io to io do if (o<>0) then begin
	 k:=dperp_offset(-nu_offs,nu[l]);
	 if node_z_alligned(l,o)=0 then begin
	    nsf[k]:=nsf[k]+1;
	    {write(l:1,' ',ltc[i,o]:1,' ');
	    if o=1 then writeln;}
	 end;
      end;
   end;
   write('s-faults: ');
   j:=0; for o:=-tfla to tfla do if (nsf[o]>0) then begin
      j:=j+1;
      write(' ',nsf[o]/2.0:6:1,'(',o:1,')');
   end;
   writeln;
   if (j=0) then writeln('none.');
end; { stacking_faults }
procedure cut_bond(Var r1,r2 : v2; Var disp:boolean);
var ox,oy,ex,ey,a,b	 : myreal;
   b1,b2,bx1,bx2,by1,by2 : boolean;
   dx,dy		 : myreal;
   o			 : integer;
begin
   x1:=r1[1]; y1:=r1[2]; x2:=r2[1]; y2:=r2[2];
   {cut the bonds at the window borders}
   {are the endpoints in or out?}
   bx1:=(x1>pswin[1]) and (x1<pswin[3]);
   bx2:=(x2>pswin[1]) and (x2<pswin[3]);
   by1:=(y1>pswin[2]) and (y1<pswin[4]);
   by2:=(y2>pswin[2]) and (y2<pswin[4]);
   b1 := bx1 and by1;
   b2 := bx2 and by2;
   dx:=x2-x1;
   dy:=y2-y1;
   {we display bond if at least one of the endpoints is inside}
   if (b1 or b2) then begin disp:=true;
      if (x1<pswin[1]) then begin
	 y1:=y2+(pswin[1]-x2)/dx*dy; x1:=pswin[1]; end;
      if (x1>pswin[3]) then begin
	 y1:=y2+(pswin[3]-x2)/dx*dy; x1:=pswin[3]; end;
      if (x2<pswin[1]) then begin
	 y2:=y1+(pswin[1]-x1)/dx*dy; x2:=pswin[1]; end;
      if (x2>pswin[3]) then begin
	 y2:=y1+(pswin[3]-x1)/dx*dy; x2:=pswin[3]; end;
      if (y1<pswin[2]) then begin
	 x1:=x2+(pswin[2]-y2)/dy*dx; y1:=pswin[2]; end;
      if (y1>pswin[4]) then begin
	 x1:=x2+(pswin[4]-y2)/dy*dx; y1:=pswin[4]; end;
      if (y2<pswin[2]) then begin
	 x2:=x1+(pswin[2]-y1)/dy*dx; y2:=pswin[2]; end;
      if (y2>pswin[4]) then begin
	 x2:=x1+(pswin[4]-y1)/dy*dx; y2:=pswin[4]; end;
   end else disp:=false;
   {the case when both pts are out is needed for follow_bond}
   if not(disp) then begin
      if (x1<pswin[1]) then x1:=pswin[1];
      if (x1>pswin[3]) then x1:=pswin[3]; 
      if (x2<pswin[1]) then x2:=pswin[1];
      if (x2>pswin[3]) then x2:=pswin[3];
      if (y1<pswin[2]) then y1:=pswin[2];
      if (y1>pswin[4]) then y1:=pswin[4];
      if (y2<pswin[2]) then y2:=pswin[2];
      if (y2>pswin[4]) then y2:=pswin[4];
   end;
   r1[1]:=x1; r1[2]:=y1; r2[1]:=x2; r2[2]:=y2;
end; { cut_bond }
procedure cut_bond_follow(Var r1,r2 : v2; Var disp:boolean);
var
   ox,oy,ex,ey,a,b	 : myreal;
   b1,b2,bx1,bx2,by1,by2 : boolean;
   x1,x2,y1,y2		 : myreal;
   xf1,xf2,yf1,yf2	 : myreal;
   dx,dy,eps		 : myreal;
   o			 : integer;
begin
   eps:=1e-7;
   x1:=r1[1]; y1:=r1[2]; x2:=r2[1]; y2:=r2[2];
   dx:=x2-x1;
   dy:=y2-y1;
   xf1:=x1; xf2:=x2; yf1:=y1; yf2:=y2; {new coordinates}
   {cut the bonds at the window borders}
   if (x1<pswin[1]) then begin
      yf1:=y2+(pswin[1]-x2)/dx*dy; xf1:=pswin[1]; end;
   if (x1>pswin[3]) then begin
      yf1:=y2+(pswin[3]-x2)/dx*dy; xf1:=pswin[3]; end;
   if (x2<pswin[1]) then begin
      yf2:=y1+(pswin[1]-x1)/dx*dy; xf2:=pswin[1]; end;
   if (x2>pswin[3]) then begin
      yf2:=y1+(pswin[3]-x1)/dx*dy; xf2:=pswin[3]; end;
   if (y1<pswin[2]) then begin
      xf1:=x2+(pswin[2]-y2)/dy*dx; yf1:=pswin[2]; end;
   if (y1>pswin[4]) then begin
      xf1:=x2+(pswin[4]-y2)/dy*dx; yf1:=pswin[4]; end;
   if (y2<pswin[2]) then begin
      xf2:=x1+(pswin[2]-y1)/dy*dx; yf2:=pswin[2]; end;
   if (y2>pswin[4]) then begin
      xf2:=x1+(pswin[4]-y1)/dy*dx; yf2:=pswin[4]; end;
   {are the endpoints in or out?}
   bx1:=(xf1>pswin[1]-eps) and (xf1<pswin[3]+eps);
   bx2:=(xf2>pswin[1]-eps) and (xf2<pswin[3]+eps);
   by1:=(yf1>pswin[2]-eps) and (yf1<pswin[4]+eps);
   by2:=(yf2>pswin[2]-eps) and (yf2<pswin[4]+eps);
   b1 := bx1 and by1;
   b2 := bx2 and by2;
   {the case when both pts are out}
   if (xf1<pswin[1]) then xf1:=pswin[1];
   if (xf1>pswin[3]) then xf1:=pswin[3]; 
   if (xf2<pswin[1]) then xf2:=pswin[1];
   if (xf2>pswin[3]) then xf2:=pswin[3];
   if (yf1<pswin[2]) then yf1:=pswin[2];
   if (yf1>pswin[4]) then yf1:=pswin[4];
   if (yf2<pswin[2]) then yf2:=pswin[2];
   if (yf2>pswin[4]) then yf2:=pswin[4];
   r1[1]:=xf1; r1[2]:=yf1; r2[1]:=xf2; r2[2]:=yf2;
   disp:=(dist(r1,r2)>eps);
end; { cut_bond }
procedure rxry(v :v2 ; Var w:v2);
var o	 : integer;
   psw_u : myreal;
   u	 : v2;
begin
   u:=v;
   if x2y then begin u[1]:=v[2]; u[2]:=-v[1]; end;
   psw_u:=psw[1];
   for o:=1 to 2 do w[o]:=over*(u[o]+rec[o])/cx*psw_u+psw[o]/2.0+pswin[o];
end; { rxry }
procedure rxry_rel(v :v2 ; Var w:v2);
var o	 : integer;
   psw_u : myreal;
   u	 : v2;
begin
   u:=v;
   if x2y then begin u[1]:=v[2]; u[2]:=v[1]; end;
   psw_u:=psw[1];
   for o:=1 to 2 do w[o]:=over*u[o]/cx*psw_u;
end; { rxry }
procedure draw_bond(u,v :v2; xw:myreal );
var display : boolean;
   o	    : integer;
begin
   cut_bond(u,v,display);
   if display then begin
      write(txfig,xw:1:2,' w ');
      write(txfig,u[1]:1:1,' ',u[2]:1:1,' ');
      writeln(txfig,v[1]:1:1,' ',v[2]:1:1,' L');
   end;
end; { draw_bond }
procedure draw_viewport(xw:myreal);
Var u,v	: array[1..4] of v2;
begin
   v[1,1]:=pswin[1]; v[1,2]:=pswin[2];
   v[2,1]:=pswin[3]; v[2,2]:=pswin[2];
   v[3,1]:=pswin[3]; v[3,2]:=pswin[4];
   v[4,1]:=pswin[1]; v[4,2]:=pswin[4];
   writeln(txfig,'gsave  %% viewport');
   writeln(txfig,v[1,1]:1:1,' ',v[1,2]:1:1,' ',v[2,1]:1:1,' ',v[2,2]:1:1,' ',v[3,1]:1:1,' ',v[3,2]:1:1,' ',v[4,1]:1:1,' ',v[4,2]:1:1,' ',' defport');
   {writeln(txfig,'newpath ',v[1,1]:1:1,' ',v[1,2]:1:1,' m ');
    writeln(txfig,v[2,1]:1:3,' ',v[2,2]:1:3,' lt ');
    writeln(txfig,v[3,1]:1:3,' ',v[3,2]:1:3,' lt ');
    writeln(txfig,'closepath clippath 1 sg fill 0 sg stroke');
   writeln(txfig,xw:1:2,' w 0.0 sg');}
   {
    writeln(txfig,v[1,1]:1:3,' ',v[1,2]:1:3,' ',v[2,1]:1:3,' ',v[2,2]:1:3,' L');
    writeln(txfig,v[2,1]:1:3,' ',v[2,2]:1:3,' ',v[3,1]:1:3,' ',v[3,2]:1:3,' L');
    writeln(txfig,v[3,1]:1:3,' ',v[3,2]:1:3,' ',v[4,1]:1:3,' ',v[4,2]:1:3,' L');
    writeln(txfig,v[4,1]:1:3,' ',v[4,2]:1:3,' ',v[1,1]:1:3,' ',v[1,2]:1:3,' L');
    }
end;
procedure draw_atom_bond(u,v :v2);
var display : boolean;
begin
   cut_bond(u,v,display);
   if display then begin
      write(txfig,u[1]:1:3,' ',u[2]:1:3,' ');
      writeln(txfig,v[1]:1:3,' ',v[2]:1:3,' L');
   end;
end; { draw_bond }
procedure cut_point(v:v2; Var disp:boolean);
begin
   disp:=((v[1]>pswin[1]) and (v[1]<pswin[3])
	  and (v[2]>pswin[2]) and (v[2]<pswin[4]));
end; { cut_points }
function rc(x : myreal):myreal;
var sz,xx : myreal;
begin
   if(x<zmin) then x:=x+1;
   sz:=(x-zmin)/(zmax-zmin);
   xx:=(min_size+sz*(max_size-min_size))*over/cx*psw[1];
   rc:=xx;
   {writeln(x:1:4,' ',sz:1:4,' ',xx:1:4,' ',max_size:1:4,' ',min_size:1:4);}
end; { rc }
function drc(x : myreal):myreal;
var sz : myreal;
begin
   sz:=(x-dzmi)/(dzmx-dzmi);
   drc:=(min_size+sz*(max_size-min_size))*over/cx*psw[1];
end; { rc }
procedure CLIP(x1,y1,x2,y2 : myreal);
begin
   write(txfig,'1.0 sg ');
   write(txfig,x1:1:2,' ',y1:1:2,' ',x2:1:2,' ',y2:1:2);
   writeln(txfig,' rectfill');
end; { CLIP }
procedure draw_unit_cell;
var
   k,j2,j3,j4,k1,k2 : integer;
   v,w,r1,r2	     : v2;
   xw		     : myreal;
   ucc		     : array[1..4] of v2;
begin
   for o:=1 to 2 do rec[o]:=0.0;
   xw:=1.5; {line thickness}
   writeln(txfig,'%%%%  Beg-Unit-Cell');
   writeln(txfig,'0 setgray 2 w');
   {writeln(txfig,'0 1 0 setrgbcolor [2 1] 0 setdash 2 w');}
   {writeln(txfig,'0.4 0.5 0.6 setrgbcolor 2.0 w');}
   ucc[1,1]:=-0.5; ucc[1,2]:=-0.5;
   ucc[2,1]:=-0.5; ucc[2,2]:= 0.5;
   ucc[3,1]:= 0.5; ucc[3,2]:= 0.5;
   ucc[4,1]:= 0.5; ucc[4,2]:=-0.5;
   for k:=1 to 4 do begin
      ortho(ucc[k],w);
      rxry(w,r1);
      write(txfig,r1[1]:1:1,' ',r1[2]:1:1,' ');
   end;
   writeln(txfig,' box stroke');
   {below draw uc throughout page}
   {
   k1:=e1; k2:=e2;
   for j1:=-k1 to k1 do
   for j2:=-k2 to k2 do
   for j3:=0 to 1 do
   for j4:=0 to 1 do if (j3<>j4) then begin
   v[1]:=j1-0.5; v[2]:=j2-0.5;
   ortho(v,w);
   rxry(w,r1);
   v[1]:=j1+j3-0.5; v[2]:=j2+j4-0.5;
   ortho(v,w);
   rxry(w,r2);
   draw_bond(r1,r2,xw);
	    end;
   }
   writeln(txfig,'0.0 sg [] 0 setdash');
   writeln(txfig,'%%%%  End-Unit-Cell');
end; { draw_unit_cell }
procedure draw_tiling(lay :integer);
var i,j,k,j1,j2,o  : integer;
   u,w,r0,r1,r2,sh : v2;
   disp,zin	   : boolean;
   rd0,rd	   : myreal;
   tn1,tn2	   : integer;
begin
   lwidth:=tbondthick*over/cx*psw[1];
   rd0:=abs(0.15*over/cx*psw[1]);
   writeln(txfig,'0.0 g');
   {writeln('draw_tiling ',advm:1,' ',lay:1,' ',tnodes1:1,' ',tnodes2:1,' ',(lay mod 2):1);}
   if(not(autoadv) and (lay mod 2 =0) and (multi>1)) then writeln(txfig,'[4 4] 0 setdash');
   for i:=tnodes1 to tnodes2 do begin
      zin:=((rt[i,3]>zmin) and (rt[i,3]<zmax));
      if inv_int then zin:=not(zin);
      {if (not(zin) and not(autoadv)) then
	 writeln(txfig,'[2 2] 0 sd') else writeln(txfig,'[] 0 sd');}
      {writeln(i:1,'  ',rt[i,3]:1:3,' ',zmin:1:3,' ',zmax:1:3);}
      for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	 {write(txfig,'[3 2] 0 sd ');}
	 if ((multi>1) and show_tbf) then begin
	    for o:=-1 to 1 do if (o<>0) then begin
	       k:=i+o*tnodes0;
	       if(k>tnodes) then k:=k-tnodes;
	       if(k<1) then k:=k+tnodes;
	       {if ((j1=0) and (j2=0)) then begin
	       writeln(tnodes:1,' ',tnodes0:1,' ',i:1,' ',k:1,' ',o:1,' ',ltm[k,-o,0]);
	       end;}
	       if (ltm[k,-o,0]=0) then begin {draw SF}
		  v[1]:=bound(rt[k,1])+j1; v[2]:=bound(rt[k,2])+j2;
		  ortho(v,r0);
		  rxry(r0,r1);
		  cut_point(r1,disp);
		  if disp then begin
		     write(txfig,r1[1]:1:1,' ',r1[2]:1:1,' ');
		     rd:=rd0/2;
		     if (o=1) then begin
			rd:=rd0/2;
			writeln(txfig,lwidth*4:1:1,' w ',rd:1:2,' C s %% SF: ',i:1,' ',o:1);
		     end else begin
			rd:=rd*0.5;
			writeln(txfig,lwidth:1:1,' w ',rd:1:2,' cro %% SF: ',i:1,' ',o:1);
		     end;
		  end;
	       end;
	    end;
	 end;
	 {writeln(txfig,' [] 0 sd');}
	 v[1]:=bound(rt[i,1])+j1; v[2]:=bound(rt[i,2])+j2;
	 ortho(v,r0);
	 rxry(r0,r1);
	 for j:=1 to 10 do begin k:=lt[i,j];
	    if (k>0) and (k>=i) then begin
	       for o:=1 to 2 do v[o]:=r0[o]+re[j,o];
	       rxry(v,r2);
	       {writeln(nu[i]:1,' ',nu[k]:1,' ',tb[nu[i],nu[k]]);}
	       if tb[nu[i],nu[k]] then begin
		  draw_bond(r1,r2,lwidth*bw[nu[i],nu[k]]);
	       end;
	    end;
	 end;
      end;
   end;
   if(not(autoadv) and (lay mod 2 =0)) then writeln(txfig,'[] 0 setdash');
end; { draw_tiling }
procedure draw_objects(beta:myreal) ;
Const
   tps	 = tdim*20;
   npstr = 800;
Type
   pstr	= string[npstr];
Var
   i,j1,j2,o,p,q,j,k,l,m       : integer;
   u,w,r0,r1,r2,sh,r0_u,r1_u   : v2;
   a			       : array[-1..1] of i10;
   ap			       : array[-1..1] of vpath;
   first,last,disp,fdisp,disp1 : boolean;
   dir			       : i10;
   np			       : array[-1..1] of integer;
   bn,ab,apb		       : array[-1..1] of boolean;
   ap0			       : vpath;
   node_passed		       : boolean;
   pat_ok,path_ok,pathok       : boolean;
   Frec,Lrec		       : pstr;
   filled		       : array[1..tps] of pstr;
   stroked		       : array[1..tps] of pstr;
   v_start,v_prev	       : v2;
   count,ang,npsF,npsL,lmax    : integer;
   stmp			       : mystring;

   procedure rlineto(Var x : pstr;y1,y2:v2);
   Var
      o	    : integer;
      xs,ys : pstr;
   begin
      ys:='';
      for o:=1 to 2 do begin
	 str(y1[o]:1:1,xs);
	 ys:=ys+xs+' ';
      end;
      for o:=1 to 2 do begin
	 str(y2[o]:1:1,xs);
	 ys:=ys+xs+' ';
      end;
      ys:=ys+' L ';      
      if(length(x)+length(ys)>=npstr)then begin
	 writeln('Insufficient array dimension npstr=',npstr:1,' in fig_.p.');
	 halt;
      end;
      x:=x+ys;
   end;
   procedure moveto(Var x : pstr;y:v2);
   Var
      o	 : integer;
      xs,ys : pstr;
   begin
      ys:='';
      for o:=1 to 2 do begin
	 str(y[o]:1:1,xs);
	 ys:=ys+xs+' ';
      end;
      ys:=ys+' m ';      
      if(length(x)+length(ys)>=npstr)then begin
	 writeln('Insufficient array dimension npstr=',npstr:1,' in fig_.p.');
	 halt;
      end;
      x:=x+ys;
   end;
   procedure lineto(Var x : pstr;y:v2);
   Var
      o	    : integer;
      xs,ys : pstr;
   begin
      ys:='';
      for o:=1 to 2 do begin
	 str(y[o]:1:1,xs);
	 ys:=ys+xs+' ';
      end;
      ys:=ys+' lt ';
      if(length(x)+length(ys)>=npstr)then begin
	 writeln('Insufficient array dimension npstr=',npstr:1,' in fig_.p.');
	 halt;
      end;
      x:=x+ys;
   end;
   procedure follow_bond(v: v2; Var nn:integer;Var pen_moved,first,last:boolean);
   var
      display	   : boolean;
      d_prev,d_now : boolean;
      o		   : integer;
      xs	   : mystring;
      vsave,vpen   : v2;
   begin
      nn:=nn+1;
      vsave:=v;
      pen_moved:=false;
      if(nn>1)then begin
	 cut_bond_follow(v_prev,v,display);
	 if(first)then begin
	    v_start:=v_prev;
	    moveto(Frec,v_start); vpen:=v_start;
	    moveto(Lrec,v_start); vpen:=v_start;
	    first:=not(first);
	 end;
	 {writeln('fb : ',d_prev,' ',d_now,' ',display);}
	 if(dist(v_prev,vpen)>0.00001)then begin
	    lineto(Frec,v_prev);
	    lineto(Lrec,v_prev);
	    vpen:=v_prev;
	    pen_moved:=true;
	 end;
	 if display then begin
	    lineto(Frec,v);
	    lineto(Lrec,v);
	    vpen:=v;
	    pen_moved:=true;
	 end;
      end;
      v_prev:=vsave;
   end; { follow_bond }
   procedure FLrecord(j	:integer);
   begin
      if (length(obj_fils[j])>1)then begin
	 npsF:=npsF+1;
	 if(npsF>tps)then begin
	    writeln('Insufficient array dim tps=',tps:1,' in fig_.p.');
	    halt;
	 end;
	 filled[npsF]:=Frec+obj_fils[j]+' %% '+obj_name[j];
      end;
      if (obj_lw[j]>0.001) then begin
	 npsL:=npsL+1;
	 if(npsL>tps)then begin
	    writeln('Insufficient array dim tps=',tps:1,' in fig_.p.');
	    halt;
	 end;
	 str(lwidth*obj_lw[j]:1:1,stmp);
	 stroked[npsL]:=Lrec+' '+stmp+' w s %% '+obj_name[j];
      end;
   end;
begin
   lmax:=0; npsF:=0; npsL:=0; {number of PS records: filled stuff, stroked stuff}
   for i:=1 to t_obj do obj_fr[i]:=0;
   lwidth:=tbondthick*over/cx*psw[1];
   if(pse) then writeln(txfig,'0.0 g % START DRAWING OBJECTS');
   for i:=tnodes1 to tnodes2 do begin
      p:=i;
      for j:=1 to n_obj do if (obj_b[0,j] or obj_b[1,j] or obj_b[2,j]) then BEGIN
	 node_passed:=true;
	 np[0]:=p; {trivial}
	 if (io3d>0) then begin
	    for o:=-io3d to io3d do if (o<>0) then begin
	       np[o]:=0;
	       bn[o]:=teq_node(ltm[p,o,0],obj_node[j,o]);
	       if (obj_node[j,o]=1) then np[o]:=ltm[p,o,0];
	    end;
	    node_passed:=(bn[-1] and bn[1]);
	 end;
	 if node_passed then for k:=1 to obj_mu[j] do if ((nu[p]=obj_nu[j,k]) or not(flavours)) then begin {loop OBJ-MULT}
	    pat_ok:=teq10(lt[p],obj_pat_list[j,k,0],a[0]);
	    if ((pat_ok) and (io3d>0)) then begin {test top/bottom pattern}
	       for o:=-io3d to io3d do if (o<>0) then begin
		  for l:=1 to 10 do begin
		     if not(teq_node(ltm[p,o,l],obj_pat_list[j,k,o,l])) then pat_ok:=false else a[o,l]:=ltm[p,o,l];
		     a[o,l]:=ltm[p,o,l];
		     if (obj_pat_list[j,k,o,l]=2) then a[o,l]:=-a[o,l];
		  end;
	       end;
	    end;
	    path_ok:=(obj_pth[j]=0);
	    if (pat_ok and not(path_ok)) then begin
	       path_ok:=fpath(p,obj_pth[j],obj_pth_def[j,k],ap0);
	       if (path_ok) then begin
		  {Note: now path MAY contain zeros!}
		  for l:=1 to obj_pth[j] do if path_ok then begin
		     {test level0}
		     path_ok:=teq_node(ap0[l],obj_pth_stat[j,0,l]);
		     if path_ok then begin
			ap[0,l]:=ap0[l];
			if (obj_pth_stat[j,0,l]=2) then ap[0,l]:=-ap[0,l];
		     end;
		     {test up/down}
		     if (path_ok and (io3d>0)) then begin
			for o:=-io3d to io3d do if ((o<>0) and path_ok) then begin
			   {ap0=0 AND path_status=2 is OK}
			   if (ap0[l]=0) then begin
			      if (obj_pth_stat[j,o,l]=2) then ap[o,l]:=0 else path_ok:=false; {vertex above nonexistent vertex is not checked}
			   end else begin
			      path_ok:=teq_node(ltm[ap0[l],o,0],obj_pth_stat[j,o,l]);
			      if path_ok then begin {keep record}
				 ap[o,l]:=ap0[l];
				 if (obj_pth_stat[j,o,l]=2) then
				    ap[o,l]:=-ap[o,l];
			      end;
			   end;
			end;
		     end; {end 3d test}
		  end; {end obj_path}
	       end;
	    end;
	    if (pat_ok and path_ok) then begin
	       obj_fr[j]:=obj_fr[j]+1;
	       {write(p:1,' ');
	       for o:=1 to 10 do write(lt[p,o]:1,','); write(' ');
	       for o:=1 to obj_pth[j] do write(ap[o]:1,' ');
	       writeln(j:1,' ',obj_pth[j]:1,' ',obj_name[j]);}
	       for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
		  v[1]:=bound(rt[i,1])+j1; v[2]:=bound(rt[i,2])+j2;
		  Frec:=''; Lrec:=''; fdisp:=false;
		  disp1:=false; {true if at least one bond cuts viewport}
		  ortho(v,r0);
		  rxry(r0,r1);
		  r0_u:=r0;
		  r1_u:=r1;
		  if obj_b[0,j] then begin {path+pattern: lines}
		     {writeln(txfig,obj_fils[j]);}
		     {follow_bond(r1,count,fdisp,first,last);}
		     for l:=1 to 10 do if (a[0,l]>0) then begin
			for o:=1 to 2 do v[o]:=r0[o]+re[l,o];
			rxry(v,r2);
			{draw_bond(r1,r2,obj_lw[j]);}
			cut_bond(r1,r2,disp);
			if disp then rlineto(Lrec,r1,r2);
		     end;
		     for l:=1 to obj_pth[j] do begin
			if (l=1) then m:=obj_pth_def[j,k,l];
			if (l>1) then m:=po(m+obj_pth_def[j,k,l]);
			for o:=1 to 2 do v[o]:=r0[o]+re[m,o];
			rxry(v,r2);
			{draw_bond(r1,r2,obj_lw[j]);}
			cut_bond(r1,r2,disp);
			if disp then rlineto(Lrec,r1,r2);
			r0:=v; r1:=r2;
		     end;
		     FLrecord(j);
		  end;
		  if obj_b[1,j] then begin {pattern as polygon}
		     count:=0; first:=true; last:=false; {when last?}
		     {follow_bond(r1,count,fdisp,first,last);}
		     m:=0;
		     for l:=1 to 10 do if (a[0,l]>0) then begin
			m:=m+1; dir[m]:=l;
			for o:=1 to 2 do v[o]:=r0[o]+re[l,o];
			rxry(v,r2);
			follow_bond(r2,count,fdisp,first,last);
			if fdisp then disp1:=true;
		     end;
		     {check whether the total angle is > Pi}
		     q:=0; for l:=2 to m do begin
			ang:=abs(dir[l]-dir[l-1]);
			if (ang>5) then q:=q+(10-ang) else q:=q+ang;
		     end;
		     {writeln(' angle : ',obj_name[j],' ',q:1);}
		     if (q>5) then begin {if yes, closepath to the 1st pt}
			for o:=1 to 2 do v[o]:=r0[o]+re[dir[1],o];
			rxry(v,r2);
			follow_bond(r2,count,fdisp,first,last);
			if fdisp then disp1:=true;
		     end else begin {otherwise, close to the current NODE}
			follow_bond(r1_u,count,fdisp,first,last);
			if fdisp then disp1:=true;
		     end;
		     follow_bond(v_start,count,fdisp,first,last);
		     if fdisp then disp1:=true;
		     if disp1 then FLrecord(j);
		  end;
		  r0:=r0_u;
		  if obj_b[2,j] then begin {path as polygon}
		     count:=0; first:=true; last:=false;
		     if(obj_start[j]=0)then begin
			follow_bond(r1,count,fdisp,first,last);
			if fdisp then disp1:=true;
		     end;
		     for l:=1 to obj_pth[j] do begin
			if (l=1) then m:=obj_pth_def[j,k,l];
			if (l>1) then m:=po(m+obj_pth_def[j,k,l]);
			for o:=1 to 2 do v[o]:=r0[o]+re[m,o];
			rxry(v,r2);
			if (l>=obj_start[j]) then begin
			   last:=(l=obj_pth[j]);
			   follow_bond(r2,count,fdisp,first,last);
			   if fdisp then disp1:=true;
			end;
			r0:=v; r1:=r2;
		     end;
		     follow_bond(v_start,count,fdisp,first,last);
		     if fdisp then disp1:=true;
		     if(disp1)then FLrecord(j);
		  end;
	       end;
	    end;
	 end;
      end; {LOOP OBJECTS}
   end; {LOOP NODES}
   lmax:=0;
   if(pse) then begin
      if(npsF>0) then writeln(txfig,'0.0 g %% objects: POLY-FILL');
      for i:=1 to npsF do begin
	 if(length(filled[i])>lmax) then lmax:=length(filled[i]);
	 writeln(txfig,filled[i]);
      end;
      if(npsL>0) then writeln(txfig,'0.0 g %% objects: LINES');
      for i:=1 to npsL do begin
	 if(length(stroked[i])>lmax) then lmax:=length(stroked[i]);
	 writeln(txfig,stroked[i]);
      end;  
      {writeln('npsF npsL MaxLength : ',npsF:1,' ',npsL:1,' ',lmax:1);}
      writeln(txfig,'%% END DRAWING OBJECTS');
   end;
   if(tinp>2) then writeln('beta=',beta:1:2);
   for i:=1 to t_obj do if (obj_fr[i]>0) then begin
      writeln(obj_fr[i]:1,' ',obj_fr[i]/tnodes:1:5,' ',obj_name[i]);
   end;
end; { draw_objects }
procedure draw_node_lab;
var i,j1,j2,o	   : integer;
   u,w,r0,r1,r2,sh : v2;
   u4		   : array[1..4] of myreal;
   disp		   : boolean;
   rd		   : myreal;
begin		   
   writeln(txfig,'0.0 g');
   fsize:=fsize0*over/cx*psw[1];
   writeln(txfig,'/Times-Roman findfont  %%% TILING');
   writeln(txfig,fsize:1:0,' scalefont');
   writeln(txfig,'setfont');
   for i:=tnodes1 to tnodes2 do begin
      for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	 v[1]:=bound(rt[i,1])+j1; v[2]:=bound(rt[i,2])+j2;
	 ortho(v,r0);
	 rxry(r0,r1);
	 cut_point(r1,disp);
	 if disp then begin
	    if (show_node_mon_mark and (fl[i]>0)) then begin
	       rd:=0.2*over/cx*psw[1];
	       writeln(txfig,' 0.0 sg 2 w ',r1[1]:1:1,' ',r1[2]:1:1,' ',rd:1:3,' C s');
	       writeln(txfig,' 1 0 0 rgb ',r1[1]:1:1,' ',r1[2]:1:1,' ',rd*0.9:1:2,' C fill');
	    end;
	    if ((show_node_mon and (fl[i]>0)) or show_node_num or show_nu or show_node_flag) then begin
	       for o:=1 to 2 do u[o]:=r1[o];
	       u4[1]:=u[1]; u4[2]:=u[2];
	       u4[3]:=0.8*fsize; u4[4]:=0.8*fsize;
	       {write(txfig,' 1.0 sg ');
		for o:=1 to 4 do write(txfig,u4[o]:1:1,' ');
		writeln(txfig,' rectfill ');
	       write(txfig,' 0.0 sg ',u[1]:1:1,' ',u[2]:1:1,' m ');}
	       o:=0;
	       if show_nu then o:=dperp_offset(-nu_offs,nu[i]);
	       if show_node_num then o:=i;
	       if show_node_mon then o:=i;
	       if show_node_flag then o:=fl[i];
	       {write(txfig,fsize*1.1:1:0,' scalefont setfont 1.0 sg');}
	       {writeln(txfig,' (',o:1,') show');}
	       {write(txfig,fsize:1:0,' scalefont setfont 0.0 sg');}
	       {writeln(txfig,' (',o:1,') show');}
	       write(txfig,' (',o:1);
	       {if(instances and (tvl[i]<0) and (tvl[i]>-ttvl))then write(txfig,tvl[i]:1,",",tvlpo[i]:1);}
	       if(instances and (tvl[i]<0) and (tvl[i]>-ttvl))then write(txfig,tvl[i]:1);
	       {if(show_node_flag) then write(txfig,fl[i]:1);}
	       write(txfig,') ',u4[1]:1:1,' ',u4[2]:1:1);
	       if(tvl[i]=-ttvl) then writeln(txfig,' strshow1') else writeln(txfig,' strshow1x');
	    end;
	 end;
      end;
   end;
end; { draw_nu }
procedure draw_atoms_node_label(u:v2;fsize:myreal;l1,l2:integer);
var x	   : myreal;
   display : boolean;
begin
   cut_point(u,display);
   if display then begin
      writeln(txfig,' (',l1:1,':',l2:1,') ',u[1]:1:1,' ',u[2]:1:1,' strshow2');
   end;
end; { draw_atoms_node_label }
procedure draw_atoms_label(u:v2;fsize:myreal;orb:integer);
var x	   : myreal;
   display : boolean;
begin
   cut_point(u,display);
   if display then begin
      writeln(txfig,' (',orb:1,') ',(u[1]):1:1,' ',(u[2]):1:1,' strshow2');
   end;
end; { draw_atoms_label }
procedure draw_atoms_z(u:v2;fsize,zc:myreal);
var x	   : myreal;
   display : boolean;
begin
   cut_point(u,display);
   if display then begin
      writeln(txfig,' (',zc:1:3,') ',u[1]:1:1,' ',u[2]:1:1,' strshow1');
   end;
end; { draw_atoms_label }
procedure draw_atoms_ene(u:v2;fsize,zc:myreal);
var x	   : myreal;
   display : boolean;
begin
   cut_point(u,display);
   if display then begin
      writeln(txfig,' (',zc:1:2,') ',u[1]:1:1,' ',u[2]:1:1,' strshow1');
   end;
end; { draw_atoms_label }
procedure draw_atom(itn,ia,itp,iorb : integer;u:v2;z:myreal;what:integer);
var
   x,xlw,rsz : myreal;
   display   : boolean;
   sh	     : v2;
   u4	     : array[1..4] of myreal;
   agr	     : mystring;
   resize,gl : myreal;
begin
   {USE atlw as a RESIZE FACTOR!!!}
   if(what=1)then begin
      xlw:=atlw[itp];
      resize:=atlw[itp];
      agr:=agrs[itp];
      if(ots[iorb])then agr:=ogrs[iorb];
   end;
   if(what=2)then begin {this is xyz}
      xlw:=atlwx[itp];
      resize:=atlwx[itp];
      agr:=agrsx[itp];
   end;
   xlw:=1.0;
   cut_point(u,display);
   if display then begin
      if (z>z2) then z2:=z;
      if (z<z1) then z1:=z;
      write(txfig,u[1]:1:1,' ',u[2]:1:1,' ');
      rsz:=rc(z); if(show_z)then rsz:=drc(n_dz[itn,ia]);
      {writeln(rsz:1:3);}
      if(itp=0) and not(show_z) then resize:=resize/3;
      write(txfig,rsz*resize:1:2,' ',atom_plot_str[itp]);
      if show_occupancy then write(txfig,n_ocr[itn,ia]:1:3,' g') else write(txfig,agr);
      writeln(txfig,' % ',itn:1,' ',ia:1,' ',itp:1,' ',iorb:1,' ',z:1:4);
      write(txfig,xlw*lwidth:1:2,' w 0.0 sg ',u[1]:1:1,' ',u[2]:1:1,' ');
      writeln(txfig,rsz*resize:1:2,atom_plot_str[itp]+'s');
      if((itp=0) and not(show_z)) then begin
	 {writeln(txfig,u[1]:1:1,' ',u[2]:1:1,' ',rsz*resize/2:1:2,' cro');}
      end;
   end;
end; { draw_atom }
procedure draw_arrow(u,v : v2);
Var
   dx,dy,x3,y3,xx	: myreal;
   width,awidth,alength	: myreal;
begin
   dx:=v[1]-u[1];
   dy:=v[2]-u[2];
   x3:=v[1]+1.07*dx;
   y3:=v[2]+1.07*dy;
   xx:=sqrt(dx*dx+dy*dy);
   if (xx>0.01) then begin
      width:=xx/20.0;
      awidth:=xx/6.0;
      alength:=xx/2.5;
      write(txfig,v[1]:1:3,' ',v[2]:1:3,' ');
      write(txfig,u[1]:1:3,' ',u[2]:1:3,' ');
      write(txfig,width*3:1:3,' ',awidth*2:1:3,' ',alength:1:3);
      writeln(txfig,' arrow 1 g fill ');
      write(txfig,v[1]:1:3,' ',v[2]:1:3,' ');
      write(txfig,u[1]:1:3,' ',u[2]:1:3,' ');
      write(txfig,width:1:3,' ',awidth:1:3,' ',alength:1:3);
      writeln(txfig,' arrow 0 g fill ');
   end;
end;
procedure draw_atom_arrow(itn,ia,itp,iorb: integer;u,v:v2;z:myreal;what:integer);
var
   x,xlw,rsz : myreal;
   display   : boolean;
   agr	     : mystring;
begin
   if(what=1)then begin
      xlw:=atlw[itp];
      agr:=agrs[itp];
      if(ots[iorb])then agr:=ogrs[iorb];
   end;
   if(what=2)then begin
      xlw:=atlwx[itp];
      agr:=agrsx[itp];
   end;
   cut_point(u,display);
   if display then begin
      {for o:=1 to 2 do write(v[o]:1:3,' ');
      for o:=1 to 2 do write(u[o]:1:3,' '); writeln;}
      if (z>z2) then z2:=z;
      if (z<z1) then z1:=z;
      if show_z then rsz:=drc(n_dz[itn,ia]) else rsz:=rc(z);
      {x:=1.0-(rho[tp]-rho_min)/(rho_max-rho_min);}
      write(txfig,u[1]:1:1,' ',u[2]:1:1,' ');
      write(txfig,rsz:1:2,' ',atom_plot_str[itp]);
      write(txfig,agr);
      writeln(txfig,' % ',itn:1,' ',ia:1,' ',itp:1,' ',iorb:1,' ',z:1:4);
      write(txfig,xlw*lwidth:1:2,' w 0.0 sg ',v[1]:1:1,' ',v[2]:1:1,' ');
      writeln(txfig,rsz:1:2,atom_plot_str[itp]+'s');
      draw_arrow(u,v);      
   end;
end; { draw_atom_arrow }
procedure draw_atom_perp(itn,ia,itp,iorb : integer;u:v2;z:myreal;what:integer);
var
   x,xlw   : myreal;
   display : boolean;
   sh	   : v2;
   u4	   : array[1..4] of myreal;
   w4	   : array[1..4] of v2;
   i,j,o   : integer;
   baps	   : array[1..2] of v2;
   xf,yf   : myreal;
   xbp	   : array[1..2] of v2;
   agr	   : mystring;
   gl	   : myreal;
begin
   if(what=1)then begin
      xlw:=atlw[itp];
      agr:=agrs[itp];
      if(ots[iorb])then begin
	 xlw:=otlw[iorb];
	 agr:=ogrs[iorb];
      end;
   end;
   if(what=2)then begin
      xlw:=atlwx[itp];
      agr:=agrsx[itp];
   end;
   if(length(agr)=0) then agr:='1 g';
   cut_point(u,display);
   if display then begin
      xf:=sqrt(5)/4.0; {need to mult by 2 for rhombus UC!}
      yf:=sqrt(5)/4.0;
      for o:=1 to 2 do xbp[1,o]:=bap[1,o]*xf;
      for o:=1 to 2 do xbp[2,o]:=bap[2,o]*yf;
      for i:=1 to 2 do rxry_rel(xbp[i],baps[i]);
      w4[1,1]:=u[1]-baps[1,1]-baps[2,1]; {--}
      w4[1,2]:=u[2]-baps[1,2]-baps[2,2];
      w4[2,1]:=u[1]+baps[1,1]-baps[2,1]; {+-}
      w4[2,2]:=u[2]+baps[1,2]-baps[2,2];
      w4[3,1]:=u[1]+baps[1,1]+baps[2,1]; {++}
      w4[3,2]:=u[2]+baps[1,2]+baps[2,2];
      w4[4,1]:=u[1]-baps[1,1]+baps[2,1]; {-+}
      w4[4,2]:=u[2]-baps[1,2]+baps[2,2];
      {alternative}
      {w4[1,1]:=u[1]-baps[2,1];
      w4[1,2]:=u[2]-baps[2,2];
      w4[2,1]:=u[1]+baps[1,1];
      w4[2,2]:=u[2]+baps[1,2];
      w4[3,1]:=u[1]+baps[2,1];
      w4[3,2]:=u[2]+baps[2,2];
      w4[4,1]:=u[1]-baps[1,1];
      w4[4,2]:=u[2]-baps[1,2];}
      if (z>z2) then z2:=z;
      if (z<z1) then z1:=z;
      for i:=1 to 4 do
	 for o:=1 to 2 do write(txfig,w4[i,o]:1:2,' ');
      write(txfig,' box ');
      if show_occupancy then write(txfig,n_ocr[itn,ia]:1:3,' g')
      else  if show_z then begin
	 gl:=(n_dz[itn,ia]-dzmi)/(dzmx-dzmi);
	 {gl:=(n_dz[itn,ia]+zlim)/(2*zlim);}
	 if(gl>1.0) then gl:=1.0;
	 if(gl<0) then gl:=0;
	 write(txfig,gl:1:3,' g')
      end
      else write(txfig,agr);
      writeln(txfig,' % ',itn:1,' ',ia:1,' ',itp:1,' ',iorb:1);
      {writeln(itn:1,' ',ia:1,' ',itp:1,' ',iorb:1,' ',otlw[iorb]:1:4,' ',xlw:1:4,' ',pboxstroke);}
      if (not(show_occupancy) and (xlw>0.001) and (pboxstroke)) then begin
	 write(txfig,xlw*lwidth:1:2,' w 0.0 sg ');
	 for i:=1 to 4 do
	    for o:=1 to 2 do write(txfig,w4[i,o]:1:2,' ');
	 writeln(txfig,' box s');
      end;
   end;
end; { draw_atom_perp }
procedure draw_perp_node(u:v2;z,xlw:myreal;nu:integer);
var x	   : myreal;
   display : boolean;
   sh	   : v2;
   u4	   : array[1..4] of myreal;
   w4	   : array[1..4] of v2;
   i,j,o   : integer;
   baps	   : array[1..2] of v2;
   xf,yf   : myreal;
   xbp	   : array[1..2] of v2;
begin
   cut_point(u,display);
   if display then begin
      xf:=sqrt(5)/4;
      yf:=sqrt(5)/4;
      for o:=1 to 2 do xbp[1,o]:=bap[1,o]*xf;
      for o:=1 to 2 do xbp[2,o]:=bap[2,o]*yf;
      for i:=1 to 2 do rxry_rel(xbp[i],baps[i]);
      w4[1,1]:=u[1]-baps[1,1]-baps[2,1];
      w4[1,2]:=u[2]-baps[1,2]-baps[2,2];
      w4[2,1]:=u[1]+baps[1,1]-baps[2,1];
      w4[2,2]:=u[2]+baps[1,2]-baps[2,2];
      w4[3,1]:=u[1]+baps[1,1]+baps[2,1];
      w4[3,2]:=u[2]+baps[1,2]+baps[2,2];
      w4[4,1]:=u[1]-baps[1,1]+baps[2,1];
      w4[4,2]:=u[2]-baps[1,2]+baps[2,2];
      if z>z2 then z2:=z;
      if z<z1 then z1:=z;
      for i:=1 to 4 do
	 for o:=1 to 2 do write(txfig,w4[i,o]:1:2,' ');
      writeln(txfig,' box s');
   end;
end; { draw_perp_node }
procedure draw_deco(plotlabel:boolean);
var
   i,j,j1,j2,o,nio,na,l	: integer;
   v,va,w,r0		: v2;
   zin,show_this	: boolean;
   zco			: myreal;
   w3			: v3;
begin
   writeln(txfig,'%%% DRAW_DECO: plotlabel ',plotlabel);
   lwidth:=tbondthick*over/cx*psw[1];
   writeln(txfig,'0.5 w');
   if (plotlabel) then
   begin
      fsize:=fsize0*over/cx*psw[1];
      writeln(txfig,'/Times-Roman findfont');
      writeln(txfig,fsize:1:0,' scalefont');
      writeln(txfig,'setfont');
   end;
{writeln('DRAW_ATOM : ',tnodes1:1,' ',tnodes2:1);}
   for i:=tnodes1 to tnodes2 do begin
      for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	 v[1]:=bound(rt[i,1])+j1; v[2]:=bound(rt[i,2])+j2;
	 ortho(v,r0);
	 nio:=0;
	 for j:=1 to nan[i] do begin
	    if (ats[n_a[i,j,1]]) then begin
	       if ((j>1) and (n_a[i,j-1,2]<>n_a[i,j,2])) then nio:=0;
	       nio:=nio+1;
	       {writeln(j:1,' ',n_a[i,j-1,2]:1,' ',n_a[i,j,2]:1,' ',nio:1);}
	       zco:=bound(n_r[i,j,3]+rt[i,3]);
	       zin:=((zco>zmin) and (zco<zmax));
	       if(not(zin)) then begin
		  zco:=zco+1;
		  zin:=((zco>zmin) and (zco<zmax));
		  if(not(zin)) then begin
		     zco:=zco-2;
		     zin:=((zco>zmin) and (zco<zmax));
		  end;
	       end;
	       {writeln(i:1,' ',j:1,' ',zco:1:5,' ',zin);}
	       if inv_int then zin:=not(zin);
	       if zin then begin
		  for o:=1 to 2 do v[o]:=n_r[i,j,o];
		  ortho(v,w);
		  for o:=1 to 2 do w[o]:=w[o]+r0[o];
		  rxry(w,v);
		  if (show_arrow and pse) then begin
		     para2xyz(n_para[i,j],n_a[i,j,1],1,0,w3);
		     for o:=1 to 2 do va[o]:=w3[o];
		     ortho(va,w);
		     for o:=1 to 2 do w[o]:=w[o]+r0[o];
		     rxry(w,va);
		  end;
		  if (plotlabel and pse) then begin
		     if show_atom_node_labels then begin
			draw_atoms_node_label(v,fsize,i,j);
		     end;
		     o:=n_a[i,j,2];
		     if show_atom_labels then o:=n_a[i,j,3];
		     if show_z_alab then draw_atoms_z(v,fsize,zco);
		     if (show_orb_labels and oldraw[i,j]) then draw_atoms_label(v,fsize,o);
		     if (show_atom_labels) then draw_atoms_label(v,fsize,o)
		  end else  
		  begin
		     o:=n_a[i,j,1];
		     if (show_xyz=1) then begin
			if (n_a[i,j,3]>nsites) then begin
			   writeln('Insufficient dimension NSITES in fig_.p.');
			   halt;
			end;
		     end;
		     if show_atoms then begin
			if ((autoadv) and (multi>1)) then begin
			   zco:=n_r[i,j,3]*multi;
			end;
			{note: we dont draw vacancy in arrow mode and "lw" is really linewidth}
			if (show_arrow and (n_a[i,j,1]>0)) then draw_atom_arrow(i,j,o,n_a[i,j,2],v,va,zco,1) else draw_atom(i,j,o,n_a[i,j,2],v,zco,1);
		     end;
		  end;
	       end;
	    end;
	 end;
      end;
   end;
end; { draw_deco }
procedure draw_deco_perp(nnu:integer);
var
   i,j,k,o,na,nac,nu_u	   : integer;
   dr,vp,vpx,vpc,vr,u,v,pc : v2;
   xx			   : v3;
   zin,nub		   : boolean;
   x5			   : v5;
   zz,xma,xmi,xw,nur	   : myreal;
   x1,y1,fsi		   : myreal;
   nu0			   : integer;
   lce			   : mystring;
begin
   na:=0;
   for o:=1 to 2 do pc[o]:=0;
   {text}
   fsi:=psw[1]/15;
   writeln(txfig,'/Times-Roman findfont');
   writeln(txfig,fsi:1:0,' scalefont');
   writeln(txfig,'setfont 0.0 sg');
   x1:=pswin[1]+0.25*fsi;
   y1:=pswin[2]+0.5*fsi;
   {write(txfig,x1:1:1,' ',y1:1:1,' m ');}
   writeln(txfig,'( nu = ',nulist[nnu]:1,' )',x1:1:1,' ',y1:1:1,' strshowperp');
   for i:=tnodes1 to tnodes2 do begin
      for j:=1 to nan[i] do begin
	 if (ats[n_a[i,j,1]]) then begin
	    for o:=1 to 5 do x5[o]:=bt[i,o]+n_r5[i,j,o];
	    v5v2p(x5,vp); {subtracts background strain}
	    na:=na+1;
	    for o:=1 to 2 do begin
	       pc[o]:=pc[o]+vp[o];
	    end;
	 end;
      end;
   end;
   for o:=1 to 2 do pc[o]:=pc[o]/na;
   nac:=0;
   nu0:=0;
   xma:=0; xmi:=10000;
   for i:=tnodes1 to tnodes2 do begin
      for j:=1 to nan[i] do begin
	 if (ats[n_a[i,j,1]]) then begin
	    for o:=1 to 5 do x5[o]:=bt[i,o]+n_r5[i,j,o];
	    v5v2p(x5,vp);
	    tra_perp(vp,vpx);
	    for o:=1 to 2 do begin
	       vp[o]:=vp[o]-pc[o];
	       if vp[o]>xma then xma:=vp[o];
	       if vp[o]<xmi then xmi:=vp[o];
	    end;
	    nur:=0; for o:=1 to 5 do nur:=nur+x5[o];
	    nu_u:=dperp_offset(-nu_offs,round(nur));
	    {writeln(nu_u:1,' ',i:1,' ',j:1,' ',vp[1]:1:5,' ',vp[2]:1:5,'   ',vpx[1]:1:1,' ',vpx[2]:1:1,'  ',rt[i,1]:1:4,' ',rt[i,2]:1:4);}
	    {if nu_u<0 then for o:=1 to 2 do vp[o]:=-vp[o];}
	    if (nu_u=0) then nu0:=nu0+1;
	    zz:=bound(n_r[i,j,3]+rt[i,3]);
	    zin:=((zz>zmin) and (zz<zmax));
	    if inv_int then zin:=not(zin);
	    nub:=((nu_u)=(nulist[nnu]));
	    if (zin and nub) then begin
	       nac:=nac+1;
	       rxry(vp,v);
	       draw_atom_perp(i,j,n_a[i,j,1],n_a[i,j,2],v,zz,1);
	       if show_arrow then begin
		  ortho(n_dxy[i,j],dr);
		  {for o:=1 to 2 do write(dr[o]:1:3,' '); writeln;}
		  for o:=1 to 2 do vpx[o]:=vp[o]+dr[o]/lpar[1]*dsca;
		  rxry(vpx,u);
		  draw_arrow(u,v)
	       end;
	    end;
	 end;
      end;
   end;
   if (nnu=1) then begin
      if (lat_centered) then lce:='centered' else lce:='primitive';
      writeln('PERP mode: center (',pc[1]:1:4,',',pc[2]:1:4,'), '+lce+' lattice:');
   end;
   writeln('nu=',nulist[nnu]:1,' (',nac:1,' sites) r_max=',xma:1:4);
   {writeln('Max-Min r_perp: ',xma:1:4,' ',xmi:1:4);}
end; { draw_deco_perp }
procedure draw_tiling_perp(nnu:integer);
var
   i,j,k,o,na,nac,nu_u : integer;
   vp,vpc,vr,u,v,pc   : v2;
   zin,nub	      : boolean;
   x5		      : v5;
   zz,xma,xmi,xw,nur  : myreal;
   x1,y1,fsi,xlw      : myreal;
   nu0		      : integer;
begin
   xlw:=0.1;
   na:=0;
   for o:=1 to 2 do pc[o]:=0;
   xma:=0; xmi:=10000;
   {text}
   fsi:=psw[1]/15;
   writeln(txfig,'/Times-Roman findfont');
   writeln(txfig,fsi:1:0,' scalefont');
   writeln(txfig,'setfont 0.0 sg');
   x1:=pswin[1]+0.5*fsi;
   y1:=pswin[2]+1.0*fsi;
   write(txfig,x1:1:1,' ',y1:1:1,' m ');
   writeln(txfig,'( tiling nu = ',nulist[nnu]:1,' ) show');
   for i:=tnodes1 to tnodes2 do begin
      for o:=1 to 5 do x5[o]:=bt[i,o];
      v5v2p(x5,vp);
      na:=na+1;
      for o:=1 to 2 do begin
	 pc[o]:=pc[o]+vp[o];
	 if (vp[o]>xma) then xma:=vp[o];
	 if (vp[o]<xmi) then xmi:=vp[o];
      end;
   end;
   for o:=1 to 2 do pc[o]:=pc[o]/na;
   {for o:=1 to 2 do write(pc[o]:1:3,' '); writeln;}
   nac:=0;
   for i:=tnodes1 to tnodes2 do begin
      for o:=1 to 5 do x5[o]:=bt[i,o];
      nur:=0; for o:=1 to 5 do nur:=nur+x5[o];
      nu_u:=dperp_offset(-nu_offs,round(nur));
      v5v2p(x5,vp);
      for o:=1 to 2 do vp[o]:=vp[o]-pc[o];
      {if nu_u<0 then for o:=1 to 2 do vp[o]:=-vp[o];}
      zz:=bound(rt[i,3]);
      zin:=((zz>zmin) and (zz<zmax));
      if inv_int then zin:=not(zin);
      nub:=((nu_u)=(nulist[nnu]));
      {write(i:1,' ',nub,' ',zin);}
      if (zin and nub) then begin
	 nac:=nac+1;
	 rxry(vp,v);
	 draw_perp_node(v,zz,xlw,nu_u);
      end;
      {writeln(' ',zz:1:4,' ',zmin:1:3,' ',zmax:1:3,' ',nac:1);}
   end;
   writeln('nu=',nulist[nnu]:1,', ',na:1,' sites, ',nac:1,' shown.');
end; { draw_tiling_perp }
procedure draw_atom_bonds;
var i,j1,j2,o,k,l,m : integer;
   t1,t2,o1,o2	    : integer;
   r3,rr	    : v3;
   u,v,w,r0,r1	    : v2;
   zin1,zin2	    : boolean;
   ar,x		    : myreal;
begin
   {we loop inefficiently for every bond first to have bond types sorted}
   for m:=1 to nab do begin
      writeln(txfig,'%%% atom_bond type ',m:1);
      writeln(txfig,agrsb[m]);
      for i:=tnodes1 to tnodes2 do
	 for l:=i to tnodes2 do
	    for j:=1 to nan[i] do
	       for k:=1 to nan[l] do begin
		  t1:=n_a[i,j,1]; t2:=n_a[l,k,1];
		  o1:=n_a[i,j,2]; o2:=n_a[l,k,2];
		  {for m:=1 to nab do}
		  if (chid[m,t1,t2] and oid[m,o1,o2]) then begin
		     for o:=1 to 3 do
			rr[o]:=bound(n_r[i,j,o]+rt[i,o]-n_r[l,k,o]-rt[l,o]);
		     ortho3(rr,r3);
		     ar:=sqrt(r3[1]*r3[1]+r3[2]*r3[2]+r3[3]*r3[3]);
		     if (abs(ar-ble[m,1])<ble[m,2]) then begin
			{writeln('TTT: ',t1:1,' ',t2:1);}
			for j1:=-e1 to e1 do
			   for j2:=-e2 to e2 do begin
			      v[1]:=bound(rt[i,1]+n_r[i,j,1])+j1;
			      v[2]:=bound(rt[i,2]+n_r[i,j,2])+j2;
			      ortho(v,w);
			      rxry(w,r0);
			      for o:=1 to 2 do v[o]:=w[o]-r3[o]/lpar[1];
			      rxry(v,r1);
			      zin1:=((bound(n_r[i,j,3]+rt[i,3])>zmin) and (bound(n_r[i,j,3]+rt[i,3])<zmax));
			      if inv_int then zin1:=not(zin1);
			      zin2:=((bound(n_r[l,k,3]+rt[l,3])>zmin) and (bound(n_r[l,k,3]+rt[l,3])<zmax));
			      if inv_int then zin2:=not(zin2);
			      {if (zin1 and zin2) then}
				 draw_atom_bond(r0,r1);
			   end;
		     end;
		  end;
	       end;
      writeln(txfig,'[] 0 sd  %%% end atom_bond ',m:1);
   end;
end; { draw_atom_bonds }

procedure load_xyz;
Var
   i,j,o,l : integer;
   zco	   : myreal;
begin
   l:=0;
   {writeln('DRAW_XYZ : ',tnodes1:1,' ',tnodes2:1);}
   for i:=1 to nsites do bxyz[i]:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
	 l:=l+1;
	 if ((i>=tnodes1) and (i<=tnodes2))then begin
	    zco:=bound(n_r[i,j,3]+rt[i,3]);
	    {zzr[l]:=n_r[i,j,3]*multi+bound(zco-rax[l,3]);}
	    zzr[l]:=bound(zco-rax[l,3]);
	    if(autoadv and (multi>1)) then begin
	       zzr[l]:=(n_r[i,j,3]+zzr[l])*multi;
	    end;
	    {writeln(i:1,' ',l:1,' ',zco:1:4,' ',rax[l,3]:1:4,' ',zzr[l]:1:4);}
	    n_dz[i,j]:=bound(zco-rax[l,3]);
	    if(n_dz[j,k]>dzmx)then begin
	       dzmx:=n_dz[j,k];
	    end;
	    if(n_dz[j,k]<dzmi)then begin
	       dzmi:=n_dz[j,k];
	    end;
	    for o:=1 to 2 do n_dxy[i,j,o]:=bound(rt[i,o]+n_r[i,j,o]-rax[l,o]);
	    if((zco>zmin) and (zco<zmax))then bxyz[l]:=1 else bxyz[l]:=0;
	    bxyz[l]:=1;
	    {writeln(l:1,' ',zco:1:3,' ',bxyz[l]:1);}
	    if inv_int then bxyz[l]:=1-bxyz[l];
	    pp[l,1]:=i;
	    pp[l,2]:=j;
	 end;
      end;
end;
procedure draw_xyz;
var
   i,j1,j2,o,l,orb : integer;
   v,w,r0	   : v2;
   zco		   : myreal;
begin
   lwidth:=tbondthick*over/cx*psw[1];
   writeln(txfig,'%%% DRAW_XYZ');
   writeln(txfig,'0.5 w');
   load_xyz;
   for i:=1 to nxyz do if (bxyz[i]=1) then begin
      zco:=bound(rax[i,3]);
      if ((autoadv) and (multi>1)) then begin
	 zco:=zzr[i];
      end;
      if (atsx[ratp[i]]) then begin
	 for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	    v[1]:=bound(rax[i,1])+j1; v[2]:=bound(rax[i,2])+j2;
	    ortho(v,w);
	    rxry(w,v);
	    o:=ratp[i];
	    j:=0; {passing info that this is flat list of atoms}
	    orb:=0; {ignore orbit info}
	    draw_atom(pp[i,1],pp[i,2],o,orb,v,zco,2);
	 end;
      end;
   end;
   if (show_atom_labels or show_z_alab_rel) then
   begin
      fsize:=fsize0*over/cx*psw[1];
      writeln(txfig,'/Times-Roman findfont');
      writeln(txfig,fsize:1:0,' scalefont');
      writeln(txfig,'setfont');
   end;
   for i:=1 to nxyz do if (bxyz[i]=1) then begin
      {writeln(i:1,' ',ratp[i]:1,' ',atsx[ratp[i]]:1);}
      if (atsx[ratp[i]]) then begin
	 for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	    v[1]:=bound(rax[i,1])+j1; v[2]:=bound(rax[i,2])+j2;
	    ortho(v,w);
	    rxry(w,v);
	    if show_atom_labels then draw_atoms_label(v,fsize,i);
	    if show_z_alab_rel then draw_atoms_z(v,fsize,rax[i,3]);
	 end;
      end;
   end;
end; { draw_xyz }
procedure draw_xyz_puck(iadv :integer);
{replacement for draw_xyz that required xyz derived from deco}
Const nnl = 20;
var
   i,j1,j2,o,l,orb,nn,kk : integer;
   v,w,r0		 : v2;
   zco,zpuck0		 : myreal;
   a1,a2,a3		 : myreal;
   z1,z2,z3,dz		 : myreal;
   size0,size1,zmi0,zma0 : myreal;
   mirsize,p0size,p1size : myreal;
   izco,sg		 : integer;
   nl			 : array[-nnl..nnl] of integer;
   rnl,pmxl,rml,zave	 : array[-nnl..nnl] of myreal;
   drawstr		 : mystring;
   plot			 : array[-nnl..nnl] of boolean;
   mirstr,pstr1,pstr2	 : array[-1..1] of mystring;
   mindz,sfa		 : myreal;
begin
   {
   NOTE: structure must be centred on a PUCKERED layer!!
   }
   mindz:=0.05*4/cparxyz;
   sfa:=lpar[1]/6.5; {normalization for symbol sizes vs tiling edge length}
   if(puckmode=1) then begin {native pucker scheme}
      {mirsize:=0.10;
      p0size:=0.08;
      p1size:=0.30;}
      mirsize:=0.05/sfa;
      p0size:=0.06/sfa;
      p1size:=0.20/sfa;
      mirstr[-1]:=' cro '; {symbol for mirror layer down}
      mirstr[1]:=' plus '; {symbol for mirror layer up}
      pstr1[-1]:=' 0.0 g '; {plotstr for puckering DOWN}
      pstr2[-1]:=' C '; {symbol for puckering DOWN}
      pstr1[0]:=' 0.6 g '; {plotstr for no puckering}
      pstr2[0]:=' C '; {symbol for no puckering}
      pstr1[1]:=' 1.0 g '; {plotstr for puckering UP}
      pstr2[1]:=' C '; {symbol for puckering UP}
   end else begin {color scheme}
      mirsize:=0.07/sfa;
      p0size:=0.08/sfa;
      p1size:=0.30/sfa;
      mirstr[-1]:=' sq '; {symbol for mirror layer down}
      mirstr[1]:=' C '; {symbol for mirror layer up}
      pstr1[-1]:='chem'; {plotstr for puckering DOWN}
      pstr1[1]:='chem'; {plotstr for puckering UP}
      pstr2[-1]:=' trid '; {symbol for puckering DOWN}
      pstr2[1]:=' tri '; {symbol for puckering UP}
   end;
   zpuck0:=zpuck;
   if not(odd(iadv)) then zpuck0:=bound(zpuck+0.5);
   for i:=-nnl to nnl do begin nl[i]:=0; plot[i]:=false; rnl[i]:=0; rml[i]:=0; pmxl[i]:=0; zave[i]:=0; end;
   plot[-1]:=true; plot[0]:=true; plot[1]:=true;
   size0:=min_size; size1:=max_size;
   z1:=100; z2:=-100;
   {writeln(zmin:1:6,' ',zmax:1:6);}
   lwidth:=tbondthick*over/cx*psw[1];
   writeln(txfig,'%%% DRAW_XYZ_PUCK');
   writeln(txfig,'0.5 w 0 g');
   {load_xyz;}
   nn:=0;
   for i:=1 to nxyz do begin
      a1:=bound(rax[i,3]-zpuck0);
      if(a1<zmin) then a1:=a1+1.0;
      if(a1>zmax) then a1:=a1-1.0;
      zco:=a1;
      if ((atsx[ratp[i]]) and (zco<zmax) and (zco > zmin)) then begin
	 nn:=nn+1;
	 if(zco<z1) then z1:=zco;
	 if(zco>z2) then z2:=zco;
      end else begin
	 {writeln(iadv:1,' ',i:1,' ',rax[i,3]:1:4,' ',zco:1:4,' ',a1:1:3,' ',nn:1);}
      end;
      if ((atsx[ratp[i]]) and (zco<zmax) and (zco > zmin)) then begin
	 zmi0:=zmin; zma0:=zmax;
	 o:=ratp[i];
	 drawstr:=agrsx[o];
	 {SPECIAL FOR PUCKERING of 8A structures}
	 {W-phase has mirror layers at 0 and +0.5}
	 izco:=round(zco*4);
	 z3:=4*zco-izco;
	 sg:=1; if(z3<0) then begin sg:=-1; z3:=-z3; end;
	 nl[izco]:=nl[izco]+1;
	 rnl[izco]:=rnl[izco]+z3*z3;
	 rml[izco]:=rml[izco]+z3*sg;
	 zave[izco]:=zave[izco]+zpuck0+zco;
	 if(z3>pmxl[izco]) then pmxl[izco]:=z3;
	 {writeln(iadv:1,' ',i:1,' ',rax[i,3]:1:4,' ',zco:1:4,' ',z3:1:4,izco:3);}
	 zmin:=0; zmax:=0.5;
	 if(izco=0) then begin
	    atom_plot_str[o]:=pstr2[sg];
	    kk:=sg; if(z3<mindz) then kk:=0;
	    if(pstr1[kk]<>'chem') then agrsx[o]:=pstr1[kk];
	    min_size:=p0size;
	    max_size:=p1size;
	 end;
	 if((izco=-1) or (izco=1)) then begin
	    atom_plot_str[o]:=mirstr[izco];
	    min_size:=mirsize;
	    max_size:=mirsize;
	 end;
	 j:=0; {passing info that this is flat list of atoms}
	 orb:=0; {ignore orbit info}
	 if(plot[izco]) then begin
	    bxyz[i]:=1;
	    for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	       v[1]:=bound(rax[i,1])+j1; v[2]:=bound(rax[i,2])+j2;
	       ortho(v,w);
	       rxry(w,v);
	       draw_atom(1,1,o,orb,v,z3,2);
	    end;
	 end;
	 {restore original settings}
	 zmin:=zmi0; zmax:=zma0; agrsx[o]:=drawstr;
      end;
   end;
   if (show_atom_labels or show_z_alab_rel) then
   begin
      fsize:=fsize0*over/cx*psw[1];
      writeln(txfig,'/Times-Roman findfont');
      writeln(txfig,fsize:1:0,' scalefont');
      writeln(txfig,'setfont');
   end;
   for i:=1 to nxyz do if (bxyz[i]=1) then begin
      {writeln(i:1,' ',ratp[i]:1,' ',atsx[ratp[i]]:1);}
      if (atsx[ratp[i]]) then begin
	 for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	    v[1]:=bound(rax[i,1])+j1; v[2]:=bound(rax[i,2])+j2;
	    ortho(v,w);
	    rxry(w,v);
	    if show_atom_labels then draw_atoms_label(v,fsize,i);
	 end;
      end;
   end;
   {writeln('draw_xyz_puck: found zmin,zmax : ',z1:1:6,',',z2:1:6);}
   writeln('cpar = ',cparxyz:1:3);
   writeln('lay nat       Z      dZ   rmsdZ    MXdZ');
   for i:=-nnl to nnl do if (nl[i]>0) then begin
      write(i:2,nl[i]:5,zave[i]/nl[i]:8:4,rml[i]/nl[i]/4*cparxyz:8:4,sqrt(rnl[i]/nl[i])/4*cparxyz:8:4,pmxl[i]/4*cparxyz:8:4);
      if(plot[i]) then writeln(' PLOTTED') else writeln;
   end;
   min_size:=size0; max_size:=size1;
end; { draw_xyz_puck }
procedure draw_xyz_raw(iadv :integer);
{replacement for draw_xyz that required xyz derived from deco}
var
   i,j1,j2,o,l,orb,nn,sg : integer;
   v,w,r0		 : v2;
   zco,delta,zshift	 : myreal;
   a1,a2,a3		 : myreal;
   z1,z2,z3,dz		 : myreal;
   zmin0,zmax0		 : myreal;
   pstr1		 : array[-1..1] of mystring;
   zerode,ene,sfa	 : myreal;
begin
   {these are needed for energy show}
   sfa:=lpar[1]/4.5; {normalization for symbol sizes vs tiling edge length}
   zerode:=0.20; {fraction of dE below which atom is neither good nor bad}
   pstr1[-1]:=' 1.0 g '; {plotstr for puckering DOWN}
   pstr1[0]:=' 0.7 g '; {plotstr for puckering DOWN}
   pstr1[1]:=' 0.0 g '; {plotstr for puckering UP}
   z1:=100; z2:=-100;
   {writeln(zmin:1:6,' ',zmax:1:6);}
   delta:=-1;
   lwidth:=tbondthick*over/cx*psw[1];
   writeln(txfig,'%%% DRAW_XYZ_RAW');
   writeln(txfig,'0.5 w');
   {load_xyz;}
   if(show_atom_elab) then begin
      max_size:=0.28/sfa; min_size:=max_size*zerode/sfa;
      atom_plot_str[13]:=" C ";
      atom_plot_str[25]:=" tri ";
      atom_plot_str[26]:=" tri ";
      atom_plot_str[27]:=" tri ";
      atom_plot_str[28]:=" trid ";
   end;
   nn:=0;
   for i:=1 to nxyz do begin
      zco:=bound(rax[i,3]);
      {if ((autoadv) and (multi>1)) then begin
      a1:=bound(rax[i,3]); if(a1<0) then a1:=a1+1.0;
      a1:=a1*multi+delta;
      if(a1-iadv<zmin) then a1:=a1+multi;
      if(a1-iadv>zmax) then a1:=a1-multi;
      zco:=a1-iadv;
      end;}
      zshift:=(iadv-1)/multi;
      a1:=bound(rax[i,3]-zshift);
      a1:=a1*multi;
      if(a1<zmin) then a1:=a1+multi;
      if(a1>zmax) then a1:=a1-multi;
      zco:=a1;
      if ((atsx[ratp[i]]) and (zco<zmax) and (zco > zmin)) then begin
	 nn:=nn+1;
	 if(zco<z1) then z1:=zco;
	 if(zco>z2) then z2:=zco;
      end;
      {writeln(iadv:1,' ',i:1,' ',rax[i,3]:1:4,' ',zco:1:4,' ',a1:1:3,' ',nn:1);}
      if ((atsx[ratp[i]]) and (zco<zmax) and (zco > zmin)) then begin
	 bxyz[i]:=1;
	 if(show_atom_elab) then begin
	    zmin0:=zmin; zmax0:=zmax;
	    zmin:=0; zmax:=1;
	    zco:=sqrt(dene[i])/sqrt(demax);
	    sg:=sgs[i];
	    if(dene[i]<demax*zerode) then begin sg:=0; zco:=demax*zerode; end;
	    agrsx[ratp[i]]:=pstr1[sg];
	    {writeln(i:3,dene[i]:8:3,sg:4,demax*zerode:8:3,zco:8:3,'  ',pstr1[sg]);}
	 end;
	 for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	    v[1]:=bound(rax[i,1])+j1; v[2]:=bound(rax[i,2])+j2;
	    ortho(v,w);
	    rxry(w,v);
	    o:=ratp[i];
	    j:=0; {passing info that this is flat list of atoms}
	    orb:=0; {ignore orbit info}
	    draw_atom(1,1,o,orb,v,zco,2);
	 end;
	 if(show_atom_elab) then begin
	    zmin:=zmin0; zmax:=zmax0;
	 end;
      end;
   end;
   writeln('iadv ',iadv:1,' : ',nn:1,' atoms.  ');
   if (show_atom_labels or show_z_alab_rel) then
   begin
      fsize:=fsize0*over/cx*psw[1];
      writeln(txfig,'/Times-Roman findfont');
      writeln(txfig,fsize:1:0,' scalefont');
      writeln(txfig,'setfont');
   end;
   for i:=1 to nxyz do if (bxyz[i]=1) then begin
      {writeln(i:1,' ',ratp[i]:1,' ',atsx[ratp[i]]:1);}
      if (atsx[ratp[i]]) then begin
	 for j1:=-e1 to e1 do for j2:=-e2 to e2 do begin
	    v[1]:=bound(rax[i,1])+j1; v[2]:=bound(rax[i,2])+j2;
	    ortho(v,w);
	    rxry(w,v);
	    if show_atom_labels then draw_atoms_label(v,fsize,i);
	    {if (show_atom_elab) then draw_atoms_z(v,fsize,sene[i]);}
	    if show_z_alab_rel then draw_atoms_z(v,fsize,rax[i,3]);
	 end;
      end;
   end;
   writeln('draw_xyz_raw: found zmin_xyz,zmax_xyz ',z1:1:6,',',z2:1:6);
   if(show_atom_elab) then writeln(' max dE=',demax:1:3);
end; { draw_xyz_raw }
procedure readword(var x : text; Var y:mystring);
Var
   a		   : char;
   done,checkspace : boolean;
   k		   : integer;
begin
   y:='';
   done:=false; checkspace:=false;
   repeat
      read(x,a);
      if(checkspace) then begin
	 if(a=' ') then done:=true;
      end;
      if(a<>' ') then begin checkspace:=true; y:=y+a; end;
   until done;
end; { readword }
procedure redef_z_interval;
var x : myreal;
begin
   x:=zmin;
   zmin:=zmax;
   zmax:=x;
   inv_int:=not(inv_int);
end; { redef_z_interval }
procedure READ_INPUTF;
var
   control		 : i2;
   x,y,x1,x2,y1,y2	 : myreal;
   txt,txyz		 : text;
   i,j,o,k,k0,ol,ndef : integer;
   k1,k2,j1,j2,nline	 : integer;
   sx			 : mystring;
   sname,psgr		 : mystring;
   match,shobj		 : boolean;
   isw			 : i3;

procedure strbegend(x : mystring;chb,che:char;var y:mystring);
var
   i,j,le : integer;
   i1,i2  : integer;
begin
   y:='';
   j:=length(x);
   i1:=0; i2:=0;
   i:=0; repeat i:=i+1;
      if (x[i]=chb) then i1:=i;
   until ((i1>0) or (i=j));
   if (i1>0) then begin
      repeat i:=i+1;
	 if (x[i]=che) then i2:=i;
      until ((i2>0) or (i=j));
   end;
   if ((i1=0) or (i2=0)) then begin
      writeln('String delimiters ',chb,che,' unpaired.');
      halt;
   end;
   y:=substr(x,i1+1,(i2-i1-1));
end;
procedure startread(Var txt : text; Var x:char; Var eofi:boolean);
begin
   eofi:=false;
   repeat
      if (eoln(txt)) then begin nline:=nline+1; readln(txt); end;
      eofi:=eof(txt);
   until not(eoln(txt)) or eofi;
   if not(eofi) then begin
      read(txt,x);
      if ((x<>'=') and (x<>'<')) then begin
	 readln(txt); eofi:=eof(txt);
	 if not(eofi) then startread(txt,x,eofi);
      end;
   end;
end; { startread }
begin
   nline:=0; shobj:=false;
   reset(txt,'INPUT/fig.inp');
   readln(txt); {comment line}
   readln(txt,psout);
   readln(txt,o); title:=(o=1);
   caption:='';
   if(title) then begin
      reset(txp,'INPUT/figcap.inp');
      readln(txp,caption);
      close(txp);
   end;
   read(txt,o); show_perp:=(o=1); {perp}
   if (show_perp) then begin
      read(txt,nnup);
      if (nnup>5) then begin
	 writeln('Max 5 discrete perp levels. Setting nnup=5.');
	 nnup:=5;
      end;
      for i:=1 to nnup do read(txt,nulist[i]);
   end;
   readln(txt);
   readln(txt,o); pse:=(o=0); {nofig (no ps export)}
   readln(txt,aspr); if(aspr<0.0001) then asprb:=false else asprb:=true;
   readln(txt,o); show_tiling:=(o=1);
   if (o=1) then begin {show_tiling}
      for j1:=-tfla to tfla do
	 for j2:=-tfla to tfla do tb[j1,j2]:=true;
      if (show_perp) then show_tiling_perp:=true; {show tiling when perp is on}
   end;
   readln(txt,o,i); show_atom_elab:=(o=1);
   ecol:=0; if(show_atom_elab)  then ecol:=i;
   read(txt,show_xyz); {show_xyz}
   for i:=1 to 110 do ea[i]:=0; {mean energies} 
   for i:=1 to 110 do nea[i]:=0; {mean energies} 
   if (show_xyz>0) then begin
      read(txt,i,y,sx); readln(txt);
      readln(txt,isw[1],isw[2],isw[3]); {proj_xyz}
      strbegend(sx,'{','}',fxyz);
      if(show_xyz=3) then begin
	 write('... reading center of a puckering layer [option -pxyz]...');
	 readln(txt,zpuck,puckmode);
	 writeln(' =',zpuck:1:4);
      end;
      write('Reading xyz from '+fxyz+': ');
      reset(txyz,fxyz);
      for i:=1 to 3 do begin
	 for o:=1 to 3 do read(txyz,bl3[isw[i],o]); readln(txyz);
      end;
      cparxyz:=absv3(bl3[3]);
      readln(txyz,nats);
      if (nats>nsites) then begin
	 writeln('Insufficient array dim for '+fxyz+'.');
	 writeln('Increase Const nsites in common.p .');
	 halt;
      end;
      nxyz:=0;
      for i:=1 to nats do begin
	 nxyz:=nxyz+1;
	 for o:=1 to 3 do read(txyz,rax[nxyz,isw[o]]);
	 read(txyz,ratp[nxyz]);
	 if(ecol>0) then begin
	    read(txyz,o); {orbit #}
	    if(ecol>6) then begin
	       readword(txyz,fxyz);
	       {writeln(i:4,' ',fxyz);}
	    end;
	    read(txyz,sene[nxyz]);
	 end;
	 readln(txyz);
	 if (ratp[nxyz]=0) then nxyz:=nxyz-1 else begin
	    nea[ratp[nxyz]]:=nea[ratp[nxyz]]+1;
	    ea[ratp[nxyz]]:=ea[ratp[nxyz]]+sene[nxyz];
	 end;;
      end;
      for i:=1 to 110 do if nea[i]>0 then ea[i]:=ea[i]/nea[i];
      if(show_atom_elab) then begin
	 demax:=0;
	 demin:=10000;
	 for i:=1 to nxyz do begin
	    dene[i]:=sene[i]-ea[ratp[i]];
	    sgs[i]:=1; if(dene[i]<0) then begin dene[i]:=-dene[i]; sgs[i]:=-1; end;
	    if (dene[i]>demax) then demax:=dene[i];
	    if (dene[i]<demin) then demin:=dene[i];
	    {writeln(i:4,ratp[i]:4,dene:8:3,sgs[i]:3);}
	 end;
      end;
      close(txyz);
      writeln(nxyz:1,' atoms.');
   end else readln(txt); {end show_xyz_input}
   readln(txt,addobj);
   if(addobj=1)then begin
      writeln('Adding objects from INPUT/objects.fig.');
      READ_DECO(addobj); {add objects from objects.fig}
   end;
   read(txt,o); show_occupancy:=(o=1); {show_occupancy}
   if (show_occupancy) then begin
      read(txt,o);
      show_at_oc[o]:=true;
   end;
   readln(txt);
   readln(txt,o); show_tbf:=(o=1); {sf: show top/bottom mismatched tilenodes}
   readln(txt,o); show_orb_labels:=(o=1); {olab: orbit labels}
   readln(txt,o); show_atom_labels:=(o=1); {alab: atom labels}
   readln(txt,o); show_atom_node_labels:=(o=1); {nalab: atom-node labels}
   read(txt,o); show_z:=(o=1);
   if(show_z) then readln(txt,zlim) else readln(txt);
   readln(txt,o);
   show_z_alab:=(o=1);
   show_z_alab_rel:=(o=2);
   readln(txt,o); show_nu:=(o=1); {nnu: node nu}
   readln(txt,o); show_node_num:=(o=1); {nlab: node #s}
   readln(txt,o); x2y:=(o=1); {x2y: swap x<->y}
   readln(txt,advm); {adv: advance viewport every advm-th tiling layer}
   if (advm<0) then advm:=0;
   if (advm=0) then autoadv:=false;
   if (advm=1) then autoadv:=true else autoadv:=false;
   readln(txt,o); show_uc:=(o=1);
   readln(txt,tbondthick); {lw: default linewidth - tiling}
   readln(txt,bview); {bview: viewport border}
   readln(txt,min_size,max_size); {asizes: atom-circle radii}
   readln(txt,o); pboxstroke:=(o=1); {pboxstroke}
   readln(txt,zmin,zmax); {zwin: z-slice}
   zmin:=zmin+0.001; zmax:=zmax+0.001;
   if (zmin>zmax) then redef_z_interval;
   readln(txt,fsize0); {fsize: font_size; olx DISABLED}
   readln(txt,over); {zoom}
   readln(txt,nrows,ncols); {arrange: arrange viewport}
   readln(txt,epsw); {gap: gap between viewports}
   {for o:=1 to 4 do read(txt,pswin[o]); readln(txt); DISABLED}
   {for o:=1 to 4 do read(txt,psclip[o]); readln(txt); DISABLED}
   readln(txt,e1,e2); {uc: supercell}
   writeln('supercell ',e1:1,' ',e2:1);
   readln(txt,rec0[1],rec0[2]); {?? shift figure center}
   readln(txt,zrec0[1],zrec0[2]); {?? shift tiling bonds for successive layers}
   readln(txt,show0,show1); {show-multi}
   read(txt,o); if (o=1) then begin {superimpose temperature loops}
      for o:=1 to tloops do advance[o]:=false;
      read(txt,nadv);
      for o:=1 to nadv do begin
	 read(txt,i); advance[i]:=true;
      end;
   end; readln(txt); { fig_ }
   
   read(txt,o); if (o=1) then read(txt,goback); {repeat tiling #}
   readln(txt);
   {readln(txt,o); if (o=1) then begin autodouble:=true; dbl:=2; end;}
   readln(txt,o); {flag}
   if (o=1) then begin show_mon:=true; show_node_mon:=true; show_node_mon_mark:=true; end;
   {readln(txt,o); if (o=1) then begin
   show_mon:=true; show_node_mon_mark:=true;
   end;
   readln(txt,o); if (o=1) then show_tbf:=true; 
   read(txt,o); if (o=1) then read(txt,tbfr[1],tbfr[2]); readln(txt);}
   readln(txt,o); show_arrow:=(o=1); {arrow}
   readln(txt,ndef);
   show_atoms:=(ndef>0);
   {writeln('Show_atoms: ',ndef:1); }
   for j:=1 to ndef do begin
      read(txt,i,y,sx); readln(txt);
      {writeln('atom ',i:1,' ',y:1:3,' ',sx);}
      atom_plot_str[i]:=' C ';
      if (i<0) then begin
	 i:=abs(i);
	 atom_plot_str[i]:=' tri ';
      end;
      ats[i]:=true;
      strbegend(sx,'{','}',agrs[i]);
      atlw[i]:=y;
   end;
   if (show_xyz>0) then begin
      readln(txt,ndef);
      for j:=1 to ndef do begin
	 read(txt,i,y,sx); readln(txt);
	 atom_plot_str[i]:=' C ';
	 if (i<0) then begin
	    i:=abs(i);
	    atom_plot_str[i]:=' tri ';
	 end;
	 atsx[i]:=true;
	 strbegend(sx,'{','}',agrsx[i]);
	 atlwx[i]:=y;
      end;
   end;
   readln(txt,ndef); {orbit psdef macros}
   for j:=1 to ndef do begin
      read(txt,i,y,sx); readln(txt);
      if (j<=t_orb) then begin
	 ots[i]:=true;
	 strbegend(sx,'{','}',ogrs[i]);
	 otlw[i]:=y;
      end;
   end;
   readln(txt,ndef); show_tbonds:=(ndef>0); {tbonds}
   {writeln('Show_tbonds: ',ndef:1);}
   for j:=1 to ndef do begin
      read(txt,j1,j2,x);
      k1:=dperp_offset(nu_offs,j1);
      k2:=dperp_offset(nu_offs,j2);
      {writeln(j1:1,' ',j2:1,' ',x:1:4,'  ',k1:1,' ',k2:1);}
      bw[k1,k2]:=x; bw[k2,k1]:=x;
      tb[k1,k2]:=true; tb[k2,k1]:=true;
      readln(txt);
   end;
   readln(txt,ndef); show_atom_bonds:=(ndef>0); {show_atom_bonds}
   for k:=1 to ndef do begin
      nab:=nab+1;
      readln(txt,k1,k2,j1,j2,x,y);
      readln(txt,agrsb[nab]);
      if ((k1<0) and (k2<0)) then begin
	 for i:=0 to z_ele do
	    for j:=0 to z_ele do begin
	       chid[nab,i,j]:=true;
	    end;
      end else begin
	 if (k1<0) then begin
	    for i:=0 to z_ele do begin
	       chid[nab,k2,i]:=true;
	       chid[nab,i,k2]:=true;
	    end;
	 end;
	 if (k2<0) then begin
	    for i:=0 to z_ele do begin
	       chid[nab,k1,i]:=true;
	       chid[nab,i,k1]:=true;
	    end;
	 end;
	 if ((k1>0) and (k2>0)) then begin
	    chid[nab,k1,k2]:=true;
	    chid[nab,k2,k1]:=true;
	 end;
      end;
      if ((k1<0) and (k2<0)) then begin
	 for i:=0 to t_orb do
	    for j:=0 to t_orb do begin
	       oid[nab,i,j]:=true;
	    end;
      end else begin
	 if ((j1<0) and (j2>0)) then begin
	    for i:=0 to t_orb do begin
	       oid[nab,j2,i]:=true;
	       oid[nab,i,j2]:=true;
	    end;
	 end;
	 if ((j2<0) and (j1>0)) then begin
	    for i:=0 to t_orb do begin
	       oid[nab,j1,i]:=true;
	       oid[nab,i,j1]:=true;
	    end;
	 end;
	 if (j1>0) and (j2>0) then begin
	    oid[nab,k1,k2]:=true;
	    oid[nab,k2,k1]:=true;
	 end;
      end;
      ble[nab,1]:=x; ble[nab,2]:=y;
   end;
   readln(txt,ndef); show_objects:=(ndef>0); {show_objects}
   for j:=1 to ndef do begin
      readln(txt,sname);
      read(txt,ol);
      if (ol=2) then read(txt,k0);
      read(txt,x1); {obj-linewidth}
      readln(txt);
      readln(txt,psgr);
      j2:=length(sname);
      i:=0; repeat i:=i+1;
	 match:=(strim(sname)=strim(obj_name[i]));
	 {if (length(sname)<=length(obj_name[i])) then begin
	 sx:=substr(strim(obj_name[i]),1,length(sname));
	 match:=(sname=sx);
	    end else match:=false;}
	 if match then begin
	    shobj:=true;
	    obj_b[ol,i]:=true;
	    obj_fils[i]:=psgr;
	    obj_lw[i]:=x1;
	    if (ol=2) then obj_start[i]:=k0;
	 end;
      until ((i=n_obj) or match);
   end;
   readln(txt,seed); {seed for random number generator}
   randomize;
   close(txt);
   if (shobj) then begin
      write('Checking for objects: ');
      for i:=1 to n_obj do
	 for j:=0 to 2 do if obj_b[j,i] then write(obj_name[i],' ');
      writeln;
      if show_perp then begin
	 reset(txt,'INPUT/olab.def');
	 readln(txt,ol);
	 for i:=1 to ol do begin
	    readln(txt,k1,k2,j1,j2);
	    if(obj_b[0,j2] or obj_b[1,j2] or obj_b[2,j2]) then begin
	       ots[k1]:=true;
	       ogrs[k1]:=obj_fils[j2];
	       otlw[k1]:=obj_lw[j2];
	    end;
	 end;
      end;
   end;
end; { read_inputf }
procedure def_viewport(x1,y1,x2,y2 : myreal);
var maxx,maxy,cmp : myreal;
begin
   if show_perp then begin
      maxx:=540.0;
      maxy:=540.0;
      pswin[1]:=30+maxx*x1;
      pswin[3]:=30+maxx*x2;
      pswin[2]:=150+maxy*y1;
      pswin[4]:=150+maxy*y2;
      psw[1]:=(x2-x1)*maxx;
      psw[2]:=(y2-y1)*maxy;
      if (psw[2]<(cy/cx*psw[1])) then begin {UC doesnt fit}
	 cmp:=(cy/cx*psw[1])/psw[2];
	 psw[1]:=psw[1]/cmp;
	 pswin[3]:=pswin[1]+psw[1];
      end;
   end else begin
      maxx:=520.0;
      maxy:=520.0*cy/cx; {nov 03: was 720}
      if(maxy>650) then begin
	 cmp:=maxy/650;
	 maxy:=650;
	 maxx:=maxx/cmp;
      end;
      if(maxy<250) then begin
	 cmp:=250/maxy;
	 maxy:=250;
	 maxx:=maxx*cmp;
      end;
      if(asprb) then begin
	 maxy:=aspr*maxx;
      end;
      pswin[1]:=50+maxx*x1;
      pswin[3]:=50+maxx*x2;
      pswin[2]:=60+maxy*y1;
      pswin[4]:=60+maxy*y2;
      psw[1]:=(x2-x1)*maxx;
      psw[2]:=(y2-y1)*maxy;
      {UC doesnt fit}
      {if (psw[2]<(cy/cx*psw[1])) then begin
	 cmp:=(cy/cx*psw[1])/psw[2];
	 psw[1]:=psw[1]/cmp;
	 pswin[3]:=pswin[1]+psw[1];
      end;}
   end;
   psclip[1]:=pswin[1]-1; psclip[2]:=pswin[2]-1;
   psclip[3]:=pswin[3]+1; psclip[4]:=pswin[4]+1;
end; { def_viewport }
procedure pagenum(var txfig:text; i,j:integer);
begin
   {usual fsize0 is about 0.2}
   def_viewport(0,0,1,1);
   writeln(txfig,'/Times-Roman findfont');
   fsize:=fsize0/0.2*12;
   writeln(txfig,fsize:1:0,' scalefont');
   writeln(txfig,'setfont 0.0 sg');
   writeln(txfig,'(P ',i:1,'/',j:1,') ',(pswin[1]+5):1:1,' ',(pswin[2]+2):1:1,' strshow2');
end; { pagenum }
procedure pfname(var txfig : text;fn:mystring);
begin
   def_viewport(0,0,1,1);
   writeln(txfig,'/Times-Roman findfont');
   fsize:=fsize0/0.2*15;
   writeln(txfig,fsize:1:0,' scalefont');
   writeln(txfig,'setfont 0.0 sg');
   if(show_atom_elab) then begin
      writeln(txfig,'(',fn,' dEmax=',demax:1:3,'eV) ',(pswin[1]+7):1:1,' ',(pswin[4]-fsize-2):1:1,' strshow2x');
   end else begin
      writeln(txfig,'(',fn,') ',(pswin[1]+7):1:1,' ',(pswin[4]-fsize-2):1:1,' strshow2x');
   end;
end; { pfname }
procedure set_defaults;
Var
   i,j,k,o : integer;
begin
   dzmi:=1000;
   dzmx:=-1000;
   show_page:=false;
   flavours:=true;
   seed:=1;
   dbg:=0;
   dsca:=50.0;
   repre5:=false;
   over:=0.999; {ReSize factor}
   e1:=1; e2:=1; {include images of the unit cell}
   for o:=1 to 2 do rec[o]:=0.0; {recenter}
   for o:=1 to 2 do rec0[o]:=0.0; {recenter}
   for o:=1 to 2 do zrec0[o]:=0.0; {recenter}
   {ATOMS: default settings}
   show_atoms:=false;
   show_xyz:=0;
   show_tbonds:=false;
   show_tiling:=false;
   show_orb_labels:=false;
   show_objects:=false;
   show_node_num:=false;
   show_mon:=false;
   show_node_mon:=false;
   show_node_mon_mark:=false;
   show_tiling_perp:=false;
   show_nu:=false;
   autoadv:=false;
   autodouble:=false; dbl:=1;
   show_z_alab:=false;
   show_z_alab_rel:=false;
   show_tbf:=true;
   inv_int:=false;
   rho_min:=0.0;
   rho_max:=1.0;
   nnup:=0;
   nfig:=1;
   nab:=0;
   nrows:=1; ncols:=1;
   show0:=1; show1:=10000;
   tbfr[1]:=1.0; tbfr[2]:=1.0;
   goback:=-1;
   for o:=1 to tloops do advance[o]:=true;
   for o:=1 to z_ele do rho[o]:=0.2; {every atom gray}
   for o:=1 to z_ele do atlw[o]:=1.0;
   for o:=1 to z_ele do atlwx[o]:=1.0;
   for o:=1 to z_ele do ats[o]:=false; {dont show any atom}
   for o:=1 to z_ele do atsx[o]:=false; {dont show any xyz_atom}
   for o:=1 to t_orb do ots[o]:=false;
   for o:=1 to t_orb do otlw[o]:=1.0; {default linewidth for stroking orbit psmacros}
   for o:=0 to 2 do
      for j:=1 to n_obj do obj_b[j,o]:=false; {dont show any object}
   for o:=1 to n_obj do obj_lw[o]:=0.0; {dont show any object}
   for i:=-tfla to tfla do
      for j:=-tfla to tfla do begin bw[i,j]:=1.0; tb[i,j]:=false; end;
   for i:=1 to abond do begin
      for j:=0 to z_ele do
	 for k:=0 to z_ele do chid[i,j,k]:=false;
      for j:=1 to t_orb do
	 for k:=1 to t_orb do oid[i,j,k]:=false;
   end;
      
   zmin:=-0.5001; zmax:=0.5001;
   min_size:=0.04; max_size:=0.08; {units of tiling edge length}
   fsize0:=0.3; olx[1]:=-0.05; olx[2]:=-0.05; {orbit labels: font, offset}
   {fsize0:=0.05; olx[1]:=-0.03; olx[2]:=-0.03; }
   x2y:=false; {dont swap x<->y}
   {TILING: default settings}
   show_bonds:=true;
   show_atom_bonds:=false;
   show_occupancy:=false;
   for o:=0 to z_ele do show_at_oc[o]:=false;
   show_atom_labels:=false;
   show_atom_node_labels:=false;
   show_perp:=false;
   tbondthick:=0.01; {in the units of tiling edge}  
   epsw:=0.01;
   nobound:=false;
   loops:=1; multi:=1;
end;
BEGIN
   {MODES: now 0=tiling; 1=config; 3=tilings; 4=tiling-atoms}
   for i:=1 to t_obj do obj_st[i]:=0;
   BASIS;
   set_defaults;
   readln(ticoname);
   readln(tinp);
   readln(o);
   read_tvl:=(o=1);
   show_node_flag:=(o=2);
   write(ticoname,' : ');
   if (tinp=0) then writeln('input type TILING (mode ',tinp:1,')');
   if (tinp=1) then writeln('input type CONFIG (mode ',tinp:1,')');
   if (tinp=3) then writeln('input type MANY_TILINGS (mode ',tinp:1,')');
   if (tinp=4) then writeln('input type MANY_CONFIGS (mode ',tinp:1,')');
   if(read_tvl) then begin
      writeln('Instance mode ON.');
   end;
   if(show_node_flag) then begin
      writeln('Node_flag mode ON.');
   end;
   READ_DECO(0);
   READ_INPUTF; {needs obj_names}
   reset(txinp,ticoname);
   if ((tinp=3) or (tinp=4)) then begin
      read(txinp,tnodes,multi,loops);
      if (tnodes>tdim) then begin
	 writeln('Insufficient array dim for tiling (',tdim:1,'<',tnodes:1,')');
	 halt;
      end;
      tnodes0:=tnodes div multi;
      {writeln(loops:1);}
      {if (tinp=3) then for i:=1 to loops do begin
	 read(txinp,beta[i]);
	 str(beta[i]:1:1,sinfo[i]);
      end;}
      readln(txinp);
   end;
   show1:=loops;
   for i:=1 to 2 do begin
      for o:=1 to 5 do read(txinp,bai[i,o]); readln(txinp);
      i5v2(bai[i],ba[i]);
      i5v2p(bai[i],bap[i]);
      {write('bar[',i:1,'] ');
      for o:=1 to 2 do write(ba[i,o]:8:3);
      write('   bap[',i:1,'] ');
      for o:=1 to 2 do write(bap[i,o]:9:5); writeln;}
   end;
   check_lat_center;
   if ((tinp=0) or (tinp=1)) then begin
      readln(txinp,tnodes,multi);
      if (tnodes>tdim) then begin
	 writeln('Insufficient array dim for tiling (',tdim:1,'<',tnodes:1,')');
	 halt;
      end;
      tnodes0:=tnodes div multi;
   end;
   UNIT_CELL(multi); {used when plotting atom-bonds in A}
   zlim:=zlim/lpar[2];
   det:=ba[1,1]*ba[2,2]-ba[1,2]*ba[2,1];
   if ((tinp=3) or (tinp=4)) then begin
      rewrite(txyz,'xyzs.fig');
      writeln(txyz,(show1-show0+1):1,' samples');
      for i:=1 to 3 do begin
	 for o:=1 to 3 do write(txyz,bl3[i,o]:10:5); writeln(txyz);
      end;      
   end;
   ortho(rec0,rec); {shift}
   tnodes1:=1;
   tnodes2:=tnodes;
   if x2y then begin
      cy:=abs(ba[1,1])+abs(ba[2,1]);
      cx:=abs(ba[1,2])+abs(ba[2,2]);
   end else begin
      cx:=abs(ba[1,1])+abs(ba[2,1]);
      cy:=abs(ba[1,2])+abs(ba[2,2]);
   end;
   if show_perp then begin
      cx:=10; cy:=10;
   end;
   def_viewport(0.0,0.0,1.0,1.0);
   {writing out PS file}
   get_psdefs(npsdefs); {PS dictionary}
   rewrite(txfig,psout);
   writeln(txfig,'%!PS-Adobe-2.0');
   write(txfig,'%%BoundingBox: ');
   for o:=1 to 4 do write(txfig,psclip[o]:1:0,' '); writeln(txfig);
   for j:=1 to npsdefs do writeln(txfig,psdef[j]); {PS dictionary}
   writeln(txfig,'%%Page: 1');
   pages:=0; irow:=1; icol:=1; ii:=0; newpage:=false;
   if (advm=0) then advm:=multi; {superimpose all0}
   {for i:=-tfla to tfla do
   for j:=-tfla to tfla do if tb[i,j] then writeln(i:1,' ',j:1,' ',tb[i,j]);
   writeln('----------------------------------');}
   repeat ii:=ii+1;
      if (ii=show0) then first:=true else first:=false;
      if not(goback=ii-1) then begin
	 if (tinp=4) then begin
	    readln(txinp,sinfo[ii]);
	    {writeln(ii:2,' ',sinfo[ii]);}
	 end;
	 for j:=-tfla to tfla do nusta[j]:=0;
	 for j:=1 to tnodes do tvl[j]:=-ttvl;
	 for j:=1 to tnodes do begin
	    for o:=1 to 5 do read(txinp,bt[j,o]);
	    if(read_tvl) then begin
	       read(txinp,tvl[j],tvlpo[j]);
	    end;
	    if show_mon then readln(txinp,fl[j]) else readln(txinp);
	    nu[j]:=nui(bt[j]); {nui01 for "01" mode}
	    nusta[nu[j]]:=nusta[nu[j]]+1;
	    {nu[j]:=nuinp[j];}
	    {writeln(j:1,'  ',nu[j]:1);}
	    if not(nobound) then bound5(bt[j],rt[j]);
	    k:=((j-1) div tnodes0);
	    rt[j,3]:=k/multi;
	    {writeln(j:1,' ',k:1,' ',rt[j,3]:1:5);}
	 end;
	 tlevels:=0; for j:=-tfla to tfla do if (nusta[j]>0) then tlevels:=tlevels+1;
	 if(tlevels=2) then sigmaflip:=true else sigmaflip:=false;
	 TABLE_INIT; {init connectivity table}
	 if (multi>1) then STACKING_FAULTS;
	 if((instances) and not(read_tvl)) then begin
	    for i:=1 to tnodes do new_pat_instance(i,0);
	 end;
	 for i:=1 to tnodes do NEW_NODE_OBJECTS(i,1);
	 if((tinp=0) and repre5)then begin
	    nat:=0; for i:=1 to tnodes do for j:=1 to nan[i] do if (n_a[i,j,1]>0) then nat:=nat+1;
	    rewrite(txt,'atoms2nodes.5d');
	    for i:=1 to 2 do begin
	       for o:=1 to 5 do write(txt,bai[i,o]:1,' ');
	       writeln(txt,'  bas',i:1);
	    end;
	    writeln(txt,nat:1,' ',multi:1);
	    for i:=1 to tnodes do
	       for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
		  for o:=1 to 5 do write(txt,(bt[i,o]+n_r5[i,j,o]):1:0,' ');
		  writeln(txt);
	       end;
	    close(txt);
	 end;
	 if (tinp=3) then begin
	    l:=0; for i:=1 to tnodes do l:=l+nan[i];
	    writeln(txyz,l:1,' sites');
	    export_xyz_body(txyz);
	 end;
	 if ((tinp=1) or (tinp=4)) then begin
	    l:=0;
	    import_params(txinp);
	    for j:=1 to tnodes do begin
	       readln(txinp,nan[j]);
	       for k:=1 to nan[j] do begin
		  for o:=1 to natomrec do read(txinp,n_para[j,k,o]);
		  read(txinp,n_a[j,k,1],n_a[j,k,2]);
		  para2xyz(n_para[j,k],n_a[j,k,1],1,1,n_r[j,k]);
		  para2xyz(n_para[j,k],n_a[j,k,1],0,1,w3);
		  n_dz[j,k]:=w3[3];
		  for o:=1 to 2 do n_dxy[j,k,o]:=w3[o];
		  if(n_dz[j,k]>dzmx)then begin
		     dzmx:=n_dz[j,k];
		     ozmx:=n_a[j,k,2];
		  end;
		  if(n_dz[j,k]<dzmi)then begin
		     dzmi:=n_dz[j,k];
		     ozmi:=n_a[j,k,2];
		  end;
		  if repre5 then begin
		     for o:=1 to 5 do read(txinp,n_r5[j,k,o]);
		  end;
		  if ((tinp=1) and show_occupancy) then begin
		     n_ocr[j,k]:=1.0;
		     read(txinp,n);
		     for o:=1 to n do begin
			read(txinp,i,cmp);
			if show_at_oc[i] then n_ocr[j,k]:=1.0-cmp;
		     end;
		  end;
		  readln(txinp);
		  {ptr to atom in a flat list}
		  if (n_a[j,k,1]>0) then begin {pointer for flat list of atoms WITHOUT vacancies!}
		     l:=l+1;
		     n_a[j,k,3]:=l;
		  end else n_a[j,k,3]:=0;
	       end;
	    end;
	    if (tinp=4) then begin writeln(txyz,l:1,' sites');
	       EXPORT_XYZ_BODY(txyz);
	    end;
	 end;
      end else begin
	 ii:=ii-1;
	 goback:=-1;
	 if (ii>0) then advance[ii]:=false;
      end;
      if(show_z) then writeln('dz_min(orb) dz_max(orb) [A] : ',dzmi*lpar[2]:1:4,'(',ozmi:1,') ',dzmx*lpar[2]:1:4,'(',ozmx:1,')');
      for j:=1 to t_orb do oldrawn[j]:=false;
      for j:=1 to tnodes do
	 for k:=1 to nan[j] do oldraw[j,k]:=false;
      for j:=1 to tnodes do begin
	 oind:=0;
	 for k:=1 to nan[j] do begin
	    oind:=oind+1;
	    {draw orbit label only for first occuring orbit}
	    if not(oldrawn[n_a[j,k,2]])then begin
	       oldraw[j,k]:=true;
	       if(oind=omu[n_a[j,k,2]])then oldrawn[n_a[j,k,2]]:=true;
	    end;
	    if(oind=omu[n_a[j,k,2]])then oind:=0;
	 end;
      end;
      {writeln('SHOWX ',show0:1,' ',show1:1,' ',ii:1,' ',advance[ii]);}
      if ((ii>=show0) and (ii<=show1)) then begin
	 nfig:=1;
	 if autodouble then nfig:=dbl;
	 if (advm>0) then begin
	    if (advm>multi) then advm:=multi;
	    nfig:=multi;
	    if (multi mod advm =0) then nfig:=nfig div advm else begin writeln('Note: advance-viewport=',advm:1,' is incommensurate with multi=',multi:1,'. Superimposing all.'); advm:=multi; nfig:=1; end;
	 end;
	 if show_perp then begin
	    llp:=0; nfig:=nfig*nnup;
	 end;
	 lll:=0; {counts multi for show_perp}
	 tpages:=nfig*loops div advm div (ncols*nrows);
	 if (tpages*ncols*nrows < nfig*loops) then tpages:=tpages+1;
	 for ll:=1 to nfig do begin
	    if (llp=nnup) then llp:=0;
	    llp:=llp+1;
	    if (llp=1) then lll:=lll+1;
	    if (newpage) then begin
	       if (title) then pagenum(txfig,pages,tpages);
	       writeln(txfig,'showpage %% page ',pages:1);
	       writeln(txfig,'%%PageTrailer ');
	       writeln(txfig,'%%Page: ',(pages+1):1);
	       newpage:=false;
	    end;
	    z1:=100; z2:=-100;
	    if (show_perp) then  begin
	       if (llp=1) then begin
		  tnodes1:=((lll-1)*tnodes0)+1;
		  tnodes2:=lll*tnodes0;
	       end
	    end else begin
	       tnodes1:=((ll-1)*tnodes0*advm)+1;
	       tnodes2:=ll*tnodes0*advm;
	    end;
	    {writeln('ll tn1 tn2 ',ll:1,' ',lll:1,' ',llp:1,' ',tnodes1:5,tnodes2:5);}
	    shift_view_t:=((ii>show0) and advance[ii-1]);
	    {note: this disables any superposition of tilings from copas.out!}
	    shift_view_m:=true;
	    {writeln('Shift_view_m Shift_view_t ',shift_view_m,' ',shift_view_t,' ',tnodes0:1,' ',tnodes1:1,' ',tnodes2:1,' ',nfig:1,' ',multi:1);}
	    if (shift_view_m or shift_view_t or first) then begin
	       x1:=(icol-1)/ncols+epsw; y1:=(nrows-irow)/nrows+epsw*0.8;
	       x2:=(icol)/ncols-epsw; y2:=(nrows-irow+1)/nrows-epsw*0.8;
	       def_viewport(x1,y1,x2,y2);
	       draw_viewport(bview);
	       {writeln('page ',pages:1,' T-loop ',ii:1,' nfig ',ll:1,' row ',irow:1,' col ',icol:1);}
	    end;
	    if autodouble then begin
	       if odd(ll) then begin
		  def_viewport(x1,y1,x2,(y1+y2)/2.0-epsw/4.0);
		  draw_viewport(bview);
	       end;
	       if not(odd(ll)) then begin
		  def_viewport(x1,(y1+y2)/2.0+epsw/4.0,x2,y2);
		  draw_viewport(bview);
	       end;
	    end;
	    if not(show_perp) then begin
	       if show_objects then draw_objects(beta[ii]);
	       if show_bonds then begin
		  pswin_u:=pswin;
		  pswin[3]:=pswin[1]+tbfr[1]*psw[1];
		  pswin[4]:=pswin[2]+tbfr[2]*psw[2];
		  {if not(autoadv) then begin
		  writeln(txfig,'%%% beg tiling ',ii:1);
		  for i:=0 to multi-1 do begin
		  tnodes1:=i*tnodes0+1;
		  tnodes2:=tnodes1+tnodes0-1;
		  for o:=1 to 2 do rece[o]:=rec0[o]+zrec0[o]*i;
		  ortho(rece,rec);
		  draw_tiling;
		     end;
		  ortho(rec0,rec);
		  tnodes1:=1; tnodes2:=tnodes;
		  end else} begin
		     writeln(txfig,'%%% beg tiling(s) ',ii:1,'/',ll:1);
		     ortho(rec0,rec);
		     if(pse) then begin
			if(advm=1) then begin
			   draw_tiling(ii);
			end else begin
			   for i:=1 to multi do begin
			      tnodes1:=((i-1)*tnodes0)+1;
			      tnodes2:=i*tnodes0;
			      draw_tiling(i);
			   end;
			   tnodes1:=1; tnodes2:=tnodes;
			end;
		     end;
		  end;
		  pswin:=pswin_u;
	       end;
	       if (show_atom_bonds and pse) then draw_atom_bonds;
	       if (show_atoms and pse) then begin
		  {decorate by objects and atoms, then plot}
		  if (tinp<3) then begin
		  end;
		  for o:=1 to nxyz do bxyz[o]:=0;
		  draw_deco(FALSE);
		  if show_atoms then begin
		     if (show_orb_labels or show_atom_labels or show_atom_node_labels or show_z_alab) then draw_deco(TRUE);
		  end;
	       end;
	       if ((show_xyz>0) and pse) then begin
		  o:=(ll-1)*advm+1;
		  if(show_xyz=1) then draw_xyz;
		  if(show_xyz=2) then draw_xyz_raw(o);
		  if(show_xyz=3) then draw_xyz_puck(o);
	       end;
	       if (show_nu or show_node_num or show_mon) then if (pse) then draw_node_lab;
	       if (show_uc and pse) then draw_unit_cell;
	    end else if (pse) then begin {show_PERP}
	       if show_atoms then begin
		  if (show_xyz>0) then load_xyz;
		  draw_deco_perp(llp);
	       end;
	       if show_tiling_perp then draw_tiling_perp(llp);
	    end;
	    if ((tinp=3) or (tinp=4)) then if (pse) then begin
	       fsize:=fsize0/cx*psw[1]*0.9;
	       if not(show_perp) then begin
		  writeln(txfig,'/Times-Roman findfont');
		  writeln(txfig,fsize:1:0,' scalefont');
		  writeln(txfig,'setfont 0.0 sg');
		  write(txfig,(pswin[1]+5.0):1:1,' ',(pswin[4]-fsize*1.05):1:1,' m ');
		  if (tinp=4) then begin
		     writeln(txfig,'(tlayer:',ll:1,'/',multi:1,' t-loop:',ii:1,'/',loops:1,' ',sinfo[ii],')',(pswin[1]+10.0):1:1,' ',(pswin[4]-fsize*1.1):1:1,' strshows');
		  end else begin
		     writeln(txfig,'(tlayer:',ll:1,'/',multi:1,' loop:',ii:1,'/',loops:1,' beta=',beta[ii]:1:2,')',(pswin[1]+5.0):1:1,' ',(pswin[4]-fsize*1.1):1:1,' strshows');
		  end;
	       end;
	       {if not(show_perp) then begin
	       write(txfig,(pswin[1]+5.0):1:1,' ',(pswin[4]-23.0):1:1,' m ');
	       writeln(txfig,' (zmin,zmax: ',z1:1:3,' ',z2:1:3,') show');
	       end;}
	    end;
	    if (autodouble) then inv_int:=not(inv_int) else begin
	       writeln(txfig,'grestore');
	       writeln(txfig,'%%% end tiling ',ii:1,'/',ll:1);
	       if (shift_view_t or shift_view_m) then if (icol=ncols) then
	       begin irow:=irow+1; icol:=0; end;
	       if (((ii=loops) and (ll=multi) and not(show_perp)) or (irow>nrows)) then begin
		  irow:=1; icol:=0; pages:=pages+1; newpage:=true;
	       end;
	       if (shift_view_m or shift_view_t) then icol:=icol+1;
	    end;
	    if(title) then pfname(txfig,caption);
	 end;
      end;
      {writeln(ii:1,' ',irow:1,' ',icol:1,' ',advance[ii]);}
   until (ii=loops);
   if ((tpages>pages) or newpage) then begin {finish last page}
      if (tpages>pages) then pages:=pages+1;
      if ((pages>1) and (title)) then pagenum(txfig,tpages,tpages);
   end;
   {writeln(txfig,'showpage %% page ',pages:1);}
   close(txinp);
   writeln(txfig,'showpage %% END');
   close(txfig);
   if(pse) then writeln(pages:1,' page(s) of output -> '+psout);
   if deco then begin
      if ((tinp=0) or (tinp=1)) then export_xyz('xyz.fig');
   end;
   if ((tinp=3) or (tinp=4)) then close(txyz);
   export_obj_st('objsta.fig');
END.
